// ****************** ANALYZE TAMIL NADU DATA ****************** //

/// author: alejandro ganimian (with support from sharnic djaker)

/// date: july 18, 2023

/// tables

	/// main-text tables
	
		* 1. summary statistics at baseline and randomization balance
		* 2. follow-up rate in endline assessments
		* 3. impact on attendance and punctuality from unannounced observations
		* 4. impact on overall time allocation from unannounced observations
		* 5. impact on endline assessments (standardized scores)
		* 6. impact on endline WAZ and HAZ scores (complete sample)
		* 7. [prepared in excel]
	
	///	appendix-a tables
	
		* a1. summary statistics on facilitators from intervention monitoring
		* a2. impact on time allocation to health and nutrition from unannounced observations
		* a3. impact on endline assessments with inverse-probability weights (standardized scores)
		* a4. heterogeneous impacts on endline assessments (standardized scores)
		* a5. impact on endline assessments by AWC attendance (standardized scores)
		* a6. impact on endline WAZ and HAZ scores (commmon sample)
		* a7. impact on endline WAZ and HAZ scores with inverse-probability weights (complete sample)
		* a8. heterogeneous impacts on endline WAZ scores
		* a9. heterogeneous impacts on endline HAZ scores
		* a10. impacts on endline WAZ scores by AWC attendance
		* a11. impacts on endline HAZ scores by AWC attendance
		* a12. regression-based association between nutrition status and learning outcomes	
		* a13. [prepared in excel]

	///	appendix-b tables	
		
		* b1. comparison of in- and out-of-sample districts
		* b2. facilitator self-reports on time allocation from intervention monitoring
		* b3. impact on endline assessments (proportion-correct scores)
		* b4. impact on endline assessments (any-correct answers)
		* b5. impact on endline assessments on the same sample and items
		* b6. impacts on WAZ scores after dropping child outliers based on residuals
		* b7. impacts on HAZ scores after dropping child outliers based on residuals
		* b8. impacts on winsorized WAZ scores
		* b9. impacts on winsorized HAZ scores

	///	appendix-c tables	
		
		* c1. raw proportion-correct scores at baseline
		
/// figures
	
	///	main-text figures
	
		* 1. benefit-cost sensitivity analysis
	
	///	appendix-a figures
	
		* a1. [not a stata graph]
		* a2. [from an another study]
		* a3. children’s age distribution by round of data collection
		* a4. AWC composite score impact estimates with outcomes imputed for missing children
		* a5. distributional treatment effects on achievement
		* a6. distributional treatment effects on WAZ scores
		* a7. distributional treatment effects on HAZ scores
		* a8. benefit-cost sensitivity analysis (log-normal discount rate)
		
		* c1. distribution of proportion-correct scores in assessments by round of data collection
 			
// ************************************************* //

version 12
clear all
set more off

global	setglobals 	"1"
global	gentables	"0"
global	gengraphs	"0"

// ****************** SET GLOBALS ****************** //

if $setglobals==1 {

/// directory globals

	global	data		/* add folder path where data sets are stored */	
	global	graphs		/* add folder path where graphs will be savedd */
	global	tables		/* add folder path where tables will be savedd */
		
}
// ****************** GENERATE TABLES ****************** //

if $gentables==1 {

///	table 1: summary statistics at baseline and randomization balance

	/// load data 1
	
		use ${data}blanaldata1, clear
	
	///	panel a
	
		balancetable treat $awc_vars using ${tables}bal_checks.tex, replace ///
		posthead("\midrule\emph{A. Anganwadi centers}\\") booktabs ///
		leftctitle(" ") ctit("Control" "Treatment" "Difference" "N") ///
		observationscolumn noobs varlabels nofoot sd(par([ ])) ///
		fe(stratum) vce(cluster digitcode)
	
	///	panel b
	
		balancetable treat $aww_vars using ${tables}bal_checks.tex, append ///
		posthead("\midrule\emph{B. Anganwadi workers}\\") ///
		booktabs leftctitle(" ") sd(par([ ])) ///
		observationscolumn noobs varlabels nonumbers nohead nofoot ///
		fe(stratum) vce(cluster digitcode)

	/// load data 2
	
		use ${data}blanaldata2, clear

	///	relabel variables
	
		lab var bas_low_waz "Underweight"
		lab var bas_low_haz "Stunted"
		
	///	panel c
	
		balancetable treat $chi_vars if in_bl==1 using ///
		${tables}bal_checks.tex, append ///
		posthead("\midrule\emph{C. Children}\\") booktabs leftctitle(" ") ///
		sd(par([ ])) observationscolumn noobs varlabels nonumbers nohead ///
		fe(stratum) vce(cluster digitcode)

///	table 2: follow-up rate in endline assessments

	/// load data 1
	
		use ${data}irtssanaldata, clear
		
	///	lab vars
	
		lab var	per_math1	"Baseline score"
		lab var	std_math1	"Baseline score"
		*ajg: i do this to place all bl coeffs on the same row
			
	/// panel a
	
		/// run analysis
	
			eststo clear
			
			qui eststo: areg fol_elawc treat, absorb(stratum) ///
			cluster(digitcode)
		
			su fol_elawc if treat==0 
			qui estadd scalar ymean=r(mean)
		
			qui eststo: areg fol_elhh treat [aw=w] ///
			if hh_sample==1|hh_sample==2, ///
			absorb(stratum) cluster(digitcode)
		
			su fol_elhh [aw=w] if treat==0 & hh_sample==1|hh_sample==2
			qui estadd scalar ymean=r(mean)
	
		///	print table
		
			qui esttab _all using ${tables}est_att.tex, replace f ///
			refcat(treat "\emph{A. Treatment}", nolabel) nolines ///
			posthead(\midrule) booktabs label se r2 drop(_cons) ///
			width(\hsize) b(%9.3f) nonotes star(* 0.1 ** 0.05 *** 0.01) ///
			stats(N ymean, fmt(%9.0f %9.3f %9.3f %9.3f) ///
			labels("N (children)" "Control mean" )) ///	
			mlabels("\shortstack{AWC\\assessment}" ///
			"\shortstack{HH\\assessment}")
	
	///	panel b
	
		/// run analysis
	
			eststo clear
		
			qui eststo: areg fol_elawc treat female bas_child_age ///
			bas_waz06 bas_haz06 std_comb1 i_female i_age i_waz i_haz ///
			i_comb, absorb(stratum) cluster(digitcode)
			
			qui test i_female i_age i_waz i_haz i_comb
			qui estadd scalar fval=r(F)
			qui estadd scalar pval=r(p)
						
			qui eststo: areg fol_elhh treat female bas_child_age bas_waz06 ///
			bas_haz06 std_comb1 i_female i_age i_waz i_haz i_comb [aw=w] ///
			if hh_sample==1|hh_sample==2, absorb(stratum) cluster(digitcode)
			
			qui test i_female i_age i_waz i_haz i_comb
			qui estadd scalar fval=r(F)
			qui estadd scalar pval=r(p)

		///	print table	
			
			qui esttab _all using ${tables}est_att.tex, append f ///
			refcat(treat "\emph{B. Treatment and baseline}", nolabel) ///
			nolines posthead(\midrule) booktabs label se r2 drop(_cons) ///
			width(\hsize) b(%9.3f) nonotes star(* 0.1 ** 0.05 *** 0.01) ///
			stats(N fval pval, fmt(%9.0f %9.3f %9.3f) ///
			labels("N (children)" "F-ratio (Interactions)" ///
			"P-value")) mlabels(none) nonumbers
	
		
///	table 3: Impact on attendance and punctuality from unannounced observations

	/// load data 1
	
		use ${data}obsanaldata1, clear
	
	/// panel a
	
		balancetable (mean if treat==0 & obs==1) (mean if treat==1 & ///
		obs==1) (diff treat if obs==1) $att_awc_vars using ///
		${tables}estitt_una_att1.tex, replace covariates($cov) ///
		vce(cluster digitcode) booktabs varlabels sd(par([ ])) nolines ///
		ctitles("Control" "Treatment" "Col. (2)-(1)") ///
		groups("AWCs" "\shortstack{Impact on\\AWCs}", pattern(1 0 1) ///
		end("\cmidrule(lr){2-3} \cmidrule(l){4-4}")) ///
		leftctitle("\emph{A. Center-level impacts}") ///
		prehead("\begin{tabular}{l*{3}c} \toprule") ///
		posthead("\midrule Share of centers that were...\\") ///
		leftcobservations("N (centers)") ///
		postfoot("\end{tabular}")
		
	///	panel b
		
		balancetable (mean if treat==0 & obs==1) (mean if treat==1 & ///
		obs==1) (diff treat if obs==1) (mean if treat==1 & obs==2) ///
		$att_aww_vars using ${tables}estitt_una_att2.tex, replace ///
		covariates($cov) vce(cluster digitcode) booktabs varlabels ///
		sd(par([ ])) nolines prehead("\begin{tabular}{l*{5}c}\midrule") ///
		posthead("\midrule Share of workers who...\\") ///
		ctitles("Control" "Treatment" "Col. (2)-(1)" "Treatment") ///
		groups("AWWs" "\shortstack{Impact on\\AWWs}" ///
		"\shortstack{ECE\\facilitators}", pattern(1 0 1 1) ///
		end("\cmidrule(lr){2-3} \cmidrule(l){4-4} \cmidrule(l){5-5}")) ///
		leftctitle("\emph{B. Worker-level impacts}") ///
		leftcobservations("N (centers)") ///
		postfoot("\bottomrule\end{tabular}")

///	tables 4: impact on time allocation from unannounced observations	

	/// load data 1
	
		use ${data}obsanaldata1, clear
	
	///	print table
	
		foreach g in ovr {	
			balancetable (mean if treat==1 & obs==2) (mean if treat==0 & ///
			obs==1) (mean if treat==1 & obs==1) (diff treat if obs==1) ///
			(mean if treat==1 & obs==3) (diff treat if obs==3) ${`g'_vars} ///
			using ${tables}estitt_una_`g'.tex, replace covariates($cov) ///
			vce(cluster digitcode) booktabs varlabels sd(par([ ])) nolines ///
			prehead("\begin{tabular}{l*{7}c}\midrule") ///
			leftctitle("Minutes per day...\\") posthead("\midrule") ///
			ctitles("Treatment" "Control" "Treatment" "Col. (3)-(2)" ///
			"Col. (1)+(3)" "Col. (5)-(2)") ///
			groups("\shortstack{ECE\\facilitators}" "AWWs" ///
			"\shortstack{Impact on\\AWWs}" ///
			"\shortstack{AWWs \&\\facilitators}" ////
			"\shortstack{Impact on\\AWCs}", pattern(1 1 0 1 1 1) ///
			end("\cmidrule(lr){2-2} \cmidrule(l){3-4} \cmidrule(l){5-5}" ///
			"\cmidrule(l){6-6}" "\cmidrule(l){7-7}")) ///
			leftcobservations("N (centers)") ///
			postfoot("\bottomrule\end{tabular}")
		}					

/// table 5: impact on endline assessments (standardized scores)

	/// load data 1
	
		use ${data}irtssanaldata, clear
		
	/// panel a
	
		/// run analysis

			eststo clear
			
			foreach s in $subjects {
				eststo est2_`s': areg stdf_`s'2 treat std_`s'1 ///
				awc_std_`s'1 $cov, absorb(stratum) cluster(digitcode)
				
				qui su std_`s'2 if treat==0
				qui estadd scalar ymean=r(mean)
				
				eststo est3_`s': areg stdf_`s'3 treat std_`s'1 ///
				awc_std_`s'1 $cov [aw=w] if (hh_sample==1| ///
				hh_sample==2), absorb(stratum) cluster(digitcode)
				
				qui su std_`s'3 [aw=w] if treat==0 & ///
				(hh_sample==1|hh_sample==2)
				qui estadd scalar ymean=r(mean)
			}
			
		///	print table
		
			qui esttab est2_* using ${tables}estitt_std2.tex, replace f ///
			refcat(treat ///
			"\emph{A. Complete sample}\\\underline{AWC assessments}", ///
			nolabel) nolines posthead(\midrule) rename(std_lang1 ///
			"std_math1" std_exec1 "std_math1" std_comb1 "std_math1") ///
			booktabs label se r2 drop(_cons std_* awc_* $cov) ///
			width(\hsize) b(%9.3f) nonotes star(* 0.1 ** 0.05 *** 0.01) ///
			mlabels("\shortstack{Math}" "\shortstack{Language}" ///
			"\shortstack{Executive\\ function}" ///
			"\shortstack{Composite\\ score}") ///
			stats(N, fmt(%9.0f) labels("N (children)"))
		
			qui esttab est3_* using ${tables}estitt_std2.tex, append f ///
			refcat(treat "\underline{HH assessments}", nolabel) ///
			nolines rename(std_lang1 "std_math1" std_exec1 "std_math1" ///
			std_comb1 "std_math1") booktabs label se r2 drop(_cons ///
			std_* awc_* $cov) width(\hsize) b(%9.3f) nonotes ///
			star(* 0.1 ** 0.05 *** 0.01) stats(N, fmt(%9.0f) ///
			labels("N (children)")) mlabels(none) nonumber
	
	///	panel b
	
		/// run analysis
	
			eststo clear
			
			foreach s in $subjects {
				qui eststo est2_`s': areg stdc_`s'2 treat std_`s'1 ///
				awc_std_`s'1 $cov if in_elawc==1 & in_elhh==1 & ///
				(hh_sample==1|hh_sample==2), ///
				absorb(stratum) cluster(digitcode)
	
				qui su stdc_`s'2 if treat==0 & in_elawc==1 & ///
				in_elhh==1 & (hh_sample==1|hh_sample==2)
				qui estadd scalar ymean=r(mean)
	
				qui eststo est3_`s': areg stdc_`s'3 treat std_`s'1 ///
				awc_std_`s'1 $cov if in_elawc==1 & in_elhh==1 & ///
				(hh_sample==1|hh_sample==2), ///
				absorb(stratum) cluster(digitcode)				

				qui su stdc_`s'3 if treat==0 & in_elawc==1 & ///
				in_elhh==1 & (hh_sample==1|hh_sample==2)
				qui estadd scalar ymean=r(mean)				
			}
			
		///	print table
		
			qui esttab est2_* using ${tables}estitt_std2.tex, append f ///
			refcat(treat ///
			"\emph{B. Common sample}\\\underline{AWC assessments}", ///
			nolabel) nolines posthead(\midrule) rename(std_lang1 ///
			"std_math1" std_exec1 "std_math1" std_comb1 "std_math1") ///
			booktabs label se drop(_cons std_* awc_* $cov) width(\hsize) ///
			b(%9.3f) nonotes star(* 0.1 ** 0.05 *** 0.01) stats(N, ///
			fmt(%9.0f) labels("N (children)")) mlabels(none) nonumber

			qui esttab est3_* using ${tables}estitt_std2.tex, append f ///
			refcat(treat "\underline{HH assessments}", nolabel) ///
			nolines rename(std_lang1 "std_math1" std_exec1 "std_math1" ///
			std_comb1 "std_math1") booktabs label se drop(_cons std_* ///
			awc_* $cov) width(\hsize) b(%9.3f) nonotes star(* 0.1 ** ///
			0.05 *** 0.01) stats(N, fmt(%9.0f) labels("N (children)")) ///
			mlabels(none) nonumber
		
	///	p-values
		
		/// prep data
				
			gen i=_n
			expand 2
			bys i: gen s1=_n==1
				
			foreach s in $subjects {
				g std_`s'=stdc_`s'2 if s1==1
				replace std_`s'=stdc_`s'3 if s1==0
						
				g bas2_std_`s'=std_`s'1*s1
				g bas3_std_`s'=std_`s'1*(1-s1)
					
				g awc2_std_`s'=awc_std_`s'1*s1
				g awc3_std_`s'=awc_std_`s'1*(1-s1)				
			}
						
			g weight=1 if s1==1
			replace weight=w if s1==0

			g aww_pass102=aww_pass10*s1
			g aww_pass103=aww_pass10*(1-s1)

			g m_pass102=m_pass10*s1
			g m_pass103=m_pass10*(1-s1)
				
			g aww_exper2=aww_exper*s1
			g aww_exper3=aww_exper*(1-s1)

			g m_exper2=m_exper*s1
			g m_exper3=m_exper*(1-s1)
				
			g treat2=treat*s1
			g treat3=treat*(1-s1)
					
			egen sample_stratum=group(s1 stratum)
		
		///	run analysis
				
			foreach s in $subjects {
				qui eststo pstd_`s': areg std_`s' s1 bas2_std_`s' ///
				bas3_std_`s' awc2_std_`s' awc3_std_`s' treat2 treat3 ///
				aww_exper2 aww_exper3 m_exper2 m_exper3 aww_pass102 ///
				aww_pass103 m_pass102 m_pass103 [aw=weight] if ///
				in_elawc==1 & in_elhh==1 & (hh_sample==1| ///
				hh_sample==2), absorb(sample_stratum) ///
				cluster(digitcode)
						
				qui test treat2=treat3
				qui estadd scalar pval=r(p)
			}	
			
		///	print table

			qui esttab pstd_* using ${tables}estitt_std2.tex, append f ///
			nolines posthead(\midrule) drop(_cons s1 bas* awc* treat* ///
			aww* m_*) booktabs label se width(\hsize) b(%9.3f) nonotes ///
			star(* 0.1 ** 0.05 *** 0.01) stats(pval, fmt(%9.3f) ///
			labels("P-value (AWC $=$ HH)")) mlabels(none) nonumber	
	
/// table 6: impact on endline WAZ/HAZ scores - full sample
		
	///	store weight globals
		
		global weight_awc	""
		global weight_hh	"[aw=w]"
	
	///	store title globals
	
		global title_waz	"WAZ score"
		global title_waz2	"Underweight\\(WAZ$<$-2)"
		global title_waz3	"Severely\\underweight\\(WAZ$<$-3)"
		global title_haz	"HAZ score"
		global title_haz2	"Stunted\\(HAZ$<$-2)"
		global title_haz3	"Severely\\stunted\\(HAZ$<$-3)"
		global title_whz	"WHZ score"
		global title_whz2	"Wasted\\(WHZ$<$-2)"
		global title_whz3	"Severely\\wasted\\(WHZ$<$-3)"

	/// load data 1
	
		use ${data}irtssanaldata, clear
	
	/// panel a
	
		///	run analysis 
		
			eststo clear	
			foreach v in waz06 low_waz very_low_waz {
				qui eststo esta_awc_`v': areg awc_`v' treat bas_waz06 ///
				$cov, absorb(stratum) cluster(digitcode)
						
				qui su awc_`v' if treat==0 & bas_waz06!=.
				qui estadd scalar ymean=r(mean)
				
				qui eststo esta_hh_`v': areg hh_`v' treat bas_waz06 ///
				$cov [aw=w] if (hh_sample==1| ///
				hh_sample==2), absorb(stratum) cluster(digitcode)
						
				qui su hh_`v' [aw=w] if treat==0 & bas_waz06!=. & ///
				(hh_sample==1|hh_sample==2)
				qui estadd scalar ymean=r(mean)	
			}	
			
		///	print table
		
			qui esttab esta_awc_* using ${tables}estitt_nut_full2.tex, ///
			replace f refcat(treat "\underline{AWC measurements}", ///
			nolabel) nolines posthead(\midrule) booktabs label se ///
			drop(_cons $cov bas*) width(\hsize) b(%9.3f) nonotes ///
			star(* 0.1 ** 0.05 *** 0.01) stats(N ymean, fmt(%9.0f %9.3f ///
			%9.3f) labels("N (children)" "Control mean")) ///
			mlabels("\shortstack{${title_waz}}" ///
			"\shortstack{${title_waz2}}" "\shortstack{${title_waz3}}") 

			qui esttab esta_hh_* using ${tables}estitt_nut_full2.tex, ///
			append f refcat(treat "\underline{HH measurements}", ///
			nolabel) nolines booktabs label se drop(_cons $cov bas*) ///
			width(\hsize) b(%9.3f) nonotes star(* 0.1 ** 0.05 *** 0.01) ///
			stats(N ymean, fmt(%9.0f %9.3f %9.3f) labels("N (children)" ///
			"Control mean")) mlabels(none) nonumbers
			
	///	panel b
			
		///	run analysis 
		
			eststo clear	
			foreach v in haz06 low_haz very_low_haz {
				qui eststo estb_awc_`v': areg awc_`v' treat bas_haz06 ///
				$cov, absorb(stratum) cluster(digitcode)
						
				qui su awc_`v' if treat==0 & bas_haz06!=.
				qui estadd scalar ymean=r(mean)
				
				qui eststo estb_hh_`v': areg hh_`v' treat bas_haz06 ///
				$cov [aw=w] if (hh_sample==1| ///
				hh_sample==2), absorb(stratum) cluster(digitcode)
						
				qui su hh_`v' [aw=w] if treat==0 & bas_haz06!=. & ///
				(hh_sample==1|hh_sample==2)
				qui estadd scalar ymean=r(mean)	
			}	
		
		///	print table
		
			qui esttab estb_awc_* using ${tables}estitt_nut_full2.tex, ///
			append f refcat(treat "\underline{AWC measurements}", ///
			nolabel) nolines posthead(\midrule) booktabs label se ///
			drop(_cons $cov bas*) width(\hsize) b(%9.3f) nonotes ///
			star(* 0.1 ** 0.05 *** 0.01) stats(N ymean, fmt(%9.0f %9.3f ///
			%9.3f) labels("N (children)" "Control mean")) ///
			mlabels("\shortstack{${title_haz}}" ///
			"\shortstack{${title_haz2}}" "\shortstack{${title_haz3}}") ///
			prehead(\toprule)

			qui esttab estb_hh_* using ${tables}estitt_nut_full2.tex, ///
			append f refcat(treat "\underline{HH measurements}", nolabel) ///
			nolines booktabs label se drop(_cons $cov bas*) width(\hsize) ///
			b(%9.3f) nonotes star(* 0.1 ** 0.05 *** 0.01) stats(N ymean, ///
			fmt(%9.0f %9.3f %9.3f) labels("N (children)" "Control mean")) ///
			mlabels(none) nonumbers

	
///	table a1: summary statistics on facilitators from intervention monitoring

	/// load data 1
	
		use ${data}imanaldata1, clear

	///	print table
	
		balancetable (mean) $intmon_vars1 using ///
		${tables}sustats_intmon1.tex, replace booktabs varlabels ///
		leftcobservations("N (centers)") nolines ///
		prehead("\begin{tabular}{l*{1}c}\toprule") sd(par([ ])) ///
		posthead("\midrule") postfoot("\bottomrule\end{tabular}")
		
///	tables a2: impact on time allocation from unannounced observations	

	/// load data 1
	
		use ${data}obsanaldata1, clear
	
	///	print table
	
		foreach g in hnn {	
			balancetable (mean if treat==1 & obs==2) (mean if treat==0 & ///
			obs==1) (mean if treat==1 & obs==1) (diff treat if obs==1) ///
			(mean if treat==1 & obs==3) (diff treat if obs==3) ${`g'_vars} ///
			using ${tables}estitt_una_`g'.tex, replace covariates($cov) ///
			vce(cluster digitcode) booktabs varlabels sd(par([ ])) nolines ///
			prehead("\begin{tabular}{l*{7}c}\midrule") ///
			leftctitle("Minutes per day...\\") posthead("\midrule") ///
			ctitles("Treatment" "Control" "Treatment" "Col. (3)-(2)" ///
			"Col. (1)+(3)" "Col. (5)-(2)") ///
			groups("\shortstack{ECE\\facilitators}" "AWWs" ///
			"\shortstack{Impact on\\AWWs}" ///
			"\shortstack{AWWs \&\\facilitators}" ////
			"\shortstack{Impact on\\AWCs}", pattern(1 1 0 1 1 1) ///
			end("\cmidrule(lr){2-2} \cmidrule(l){3-4} \cmidrule(l){5-5}" ///
			"\cmidrule(l){6-6}" "\cmidrule(l){7-7}")) ///
			leftcobservations("N (centers)") ///
			postfoot("\bottomrule\end{tabular}")
		}

/// table a3: ipw-weighted impact on endline assessments (standardized scores)

	/// load data 1
	
		use ${data}irtssanaldata, clear
		
	/// est prob of followup
	
		/// awc followup
			
			gen inv_prob2=.
			
			foreach t in 0 1 {
				probit fol_elawc female std_lang1 std_math1 std_exec1 ///
				i.stratum if treat==`t', robust
				
				predict prob2`t' if treat==`t'				
				
				replace inv_prob2=1/prob2`t' if treat==`t'
			}
		
		///	hh followup
			
			gen inv_prob3=.
				
			foreach t in 0 1 {	
				probit fol_elhh female std_lang1 std_math1 std_exec1 ///
				i.stratum if treat==`t', robust
			
				predict prob3`t' if treat==`t'
			
				replace inv_prob3=1/prob3`t' if treat==`t'
			}
			
			g combwt=w*inv_prob3
		
	/// panel a
	
		/// run analysis

			eststo clear
			
			foreach s in $subjects {
				eststo est2_`s': areg stdf_`s'2 treat std_`s'1 ///
				awc_std_`s'1 $cov [pw=inv_prob2], absorb(stratum) ///
				cluster(digitcode)
				
				qui su std_`s'2 if treat==0
				qui estadd scalar ymean=r(mean)
				
				eststo est3_`s': areg stdf_`s'3 treat std_`s'1 ///
				awc_std_`s'1 $cov [pw=combwt] if ///
				(hh_sample==1|hh_sample==2), ///
				absorb(stratum) cluster(digitcode)
				
				qui su std_`s'3 [aw=w] if treat==0 & ///
				(hh_sample==1|hh_sample==2)
				qui estadd scalar ymean=r(mean)
			}
			
		///	print table
		
			qui esttab est2_* using ${tables}estitt_std_ipw.tex, replace f ///
			refcat(treat ///
			"\emph{A. Complete sample}\\\underline{AWC assessments} (N=1514)", ///
			nolabel) nolines posthead(\midrule) rename(std_lang1 ///
			"std_math1" std_exec1 "std_math1" std_comb1 "std_math1") ///
			booktabs label se r2 drop(_cons std_* awc_* $cov) ///
			width(\hsize) b(%9.3f) nonotes star(* 0.1 ** 0.05 *** 0.01) ///
			mlabels("\shortstack{Math}" "\shortstack{Language}" ///
			"\shortstack{Executive\\ function}" ///
			"\shortstack{Composite\\ score}") stats(,)
		
			qui esttab est3_* using ${tables}estitt_std_ipw.tex, append f ///
			refcat(treat "\underline{HH assessments} (N=2075)", nolabel) ///
			nolines rename(std_lang1 "std_math1" std_exec1 "std_math1" ///
			std_comb1 "std_math1") booktabs label se r2 drop(_cons ///
			std_* awc_* $cov) width(\hsize) b(%9.3f) nonotes ///
			star(* 0.1 ** 0.05 *** 0.01) stats(,) mlabels(none) nonumber
	
	///	panel b
	
		/// run analysis
	
			eststo clear
			
			foreach s in $subjects {
				qui eststo est2_`s': areg stdc_`s'2 treat std_`s'1 ///
				awc_std_`s'1 $cov if in_elawc==1 & in_elhh==1 & ///
				(hh_sample==1|hh_sample==2) ///
				[pw=inv_prob2], absorb(stratum) cluster(digitcode)
	
				qui su stdc_`s'2 if treat==0 & in_elawc==1 & ///
				in_elhh==1 & (hh_sample==1|hh_sample==2)
				qui estadd scalar ymean=r(mean)
	
				qui eststo est3_`s': areg stdc_`s'3 treat std_`s'1 ///
				awc_std_`s'1 $cov if in_elawc==1 & in_elhh==1 & ///
				(hh_sample==1|hh_sample==2) [pw=combwt], ///
				absorb(stratum) cluster(digitcode)				

				qui su stdc_`s'3 if treat==0 & in_elawc==1 & ///
				in_elhh==1 & (hh_sample==1|hh_sample==2)
				qui estadd scalar ymean=r(mean)				
			}
			
		///	print table
		
			qui esttab est2_* using ${tables}estitt_std_ipw.tex, append f ///
			refcat(treat ///
			"\emph{B. Common sample}\\\underline{AWC assessments} (N=791)", ///
			nolabel) nolines posthead(\midrule) rename(std_lang1 ///
			"std_math1" std_exec1 "std_math1" std_comb1 "std_math1") ///
			booktabs label se drop(_cons std_* awc_* $cov) width(\hsize) ///
			b(%9.3f) nonotes star(* 0.1 ** 0.05 *** 0.01) stats(,) ///
			mlabels(none) nonumber

			qui esttab est3_* using ${tables}estitt_std_ipw.tex, append f ///
			refcat(treat "\underline{HH assessments} (N=791)", nolabel) ///
			nolines rename(std_lang1 "std_math1" std_exec1 "std_math1" ///
			std_comb1 "std_math1") booktabs label se drop(_cons std_* ///
			awc_* $cov) width(\hsize) b(%9.3f) nonotes star(* 0.1 ** ///
			0.05 *** 0.01) stats(,) mlabels(none) nonumber
		
/// table a4: heterogeneous impact on endline asssessments
	
	/// load data 1
	
		use ${data}irtssanaldata, clear
		
	///	relabel vars
		
		label var waz	"Covariate"
		label var i_waz	"Interaction"
	
	///	panel a
	
		///	run analysis
				
			eststo clear
			
			foreach i in waz female months moth10 vacancy children {
				qui eststo comb2_`i': areg stdf_comb2 treat std_comb1 ///
				`i' i_`i' awc_std_comb1 $cov, absorb(stratum) ///
				cluster(digitcode)
				
				qui eststo comb3_`i': areg stdf_comb3 treat std_comb1 ///
				`i' i_`i' awc_std_comb1 $cov [aw=w] if ///
				(hh_sample==1|hh_sample==2), ///
				absorb(stratum) cluster(digitcode)
			}
			
			qui eststo comb2_bas: areg stdf_comb2 treat std_comb1 i_comb ///
			awc_std_comb1 $cov, absorb(stratum) cluster(digitcode)
			/*ajg: can't include these in the loop because baseline comb score
			is a control variable in the first set of analyses and it is 
			the main effect in the second set of analyses */
			
			qui eststo comb3_bas: areg stdf_comb3 treat std_comb1 i_comb ///
			awc_std_comb1 $cov [aw=w] if (hh_sample==1|hh_sample==2), ///
			absorb(stratum) cluster(digitcode)
		
		///	print table
		
			qui esttab comb2_* using ${tables}hetitt_std.tex, replace f ///
			nolines posthead("\midrule" "\multicolumn{2}{l}{" ///
			"\emph{A. Complete sample}}\\" "\multicolumn{2}{l}{" ///
			"\underline{AWC measurements}}\\\addlinespace") ///
			booktabs label se b(%9.3f) nonotes star(* 0.1 ** 0.05 *** ///
			0.01) stats(N, fmt(%9.0f) labels("N (children)")) ///
			width(\hsize) drop(_cons std_comb1 awc* $cov) rename(female ///
			"waz" months "waz" moth10 "waz" vacancy "waz" children "waz" ///
			comb "waz" i_female "i_waz" i_months "i_waz" i_moth10 "i_waz" ///
			i_vacancy "i_waz" i_children "i_waz" i_comb "i_waz") ///
			mlabels("\shortstack{Baseline\\WAZ\\score}"  ////
			"Female" "\shortstack{Age\\(months)}" ///
			"\shortstack{Mother\\education}" "\shortstack{AWW\\vacancy}" ///
			"\shortstack{Class\\size}" ///
			"\shortstack{Baseline\\composite\\score}")
	
			qui esttab comb3_* using ${tables}hetitt_std.tex, append f ///
			nolines posthead("\addlinespace \multicolumn{2}{l}{" ///
			"\underline{HH measurements}}\\\addlinespace") ///
			booktabs label se b(%9.3f) nonotes star(* 0.1 ** 0.05 *** ///
			0.01) stats(N, fmt(%9.0f) labels("N (children)")) ///
			width(\hsize) drop(_cons std_comb1 awc* $cov) rename(female ///
			"waz" months "waz" moth10 "waz" vacancy "waz" children "waz" ///
			comb "waz" i_female "i_waz" i_months "i_waz" i_moth10 "i_waz" ///
			i_vacancy "i_waz" i_children "i_waz" i_comb "i_waz") ///
			mlabels(none) nonumber
		
	///	panel b
	
		/// run analysis
		
			eststo clear
		
			foreach i in waz female months moth10 vacancy children {
				qui eststo comb2_`i': areg stdc_comb2 treat std_comb1 ///
				`i' i_`i' awc_std_comb1 $cov if in_elawc==1 & ///
				in_elhh==1 & (hh_sample==1|hh_sample==2), ///
				absorb(stratum) cluster(digitcode)
				
				qui eststo comb3_`i': areg stdc_comb3 treat std_comb1 ///
				`i' i_`i' awc_std_comb1 $cov if in_elawc==1 & ///
				in_elhh==1 & (hh_sample==1|hh_sample==2), ///
				absorb(stratum) cluster(digitcode)
			}
	
			qui eststo comb2_bas: areg stdc_comb2 treat std_comb1 i_comb ///
			awc_std_comb1 $cov if in_elawc==1 & in_elhh==1 & ///
			(hh_sample==1|hh_sample==2), absorb(stratum) ///
			cluster(digitcode)
			
			qui eststo comb3_bas: areg stdc_comb2 treat std_comb1 i_comb ///
			awc_std_comb1 $cov if in_elawc==1 & in_elhh==1 & ///
			(hh_sample==1|hh_sample==2), absorb(stratum) ///
			cluster(digitcode)

		///	print table
		
			qui esttab comb2_* using ${tables}hetitt_std.tex, append f ///
			nolines posthead("\midrule" "\multicolumn{2}{l}{" ///
			"\emph{B. Common sample}}\\" "\multicolumn{2}{l}{" ///
			"\underline{AWC measurements}}\\\addlinespace") ///
			booktabs label se b(%9.3f) nonotes star(* 0.1 ** 0.05 *** ///
			0.01) width(\hsize) drop(_cons std_comb1 awc* $cov) ///
			rename(female "waz" months "waz" moth10 "waz" vacancy "waz" ///
			children "waz" comb "waz" i_female "i_waz" i_months "i_waz" ///
			i_moth10 "i_waz" i_vacancy "i_waz" i_children "i_waz" i_comb ///
			"i_waz") noobs mlabels(none) nonumber

			qui esttab comb3_* using ${tables}hetitt_std.tex, append f ///
			nolines posthead("\addlinespace \multicolumn{2}{l}{" ///
			"\underline{HH measurements}}\\\addlinespace") ///
			booktabs label se ///
			b(%9.3f) nonotes star(* 0.1 ** 0.05 *** 0.01) ///
			stats(N, fmt(%9.0f) labels("N (children)")) width(\hsize) ///
			drop(_cons std_comb1 awc* $cov) rename(female "waz" months ///
			"waz" moth10 "waz" vacancy "waz" children "waz" comb "waz" ///
			i_female "i_waz" i_months "i_waz" i_moth10 "i_waz" i_vacancy ///
			"i_waz" i_children "i_waz" i_comb "i_waz") mlabels(none) nonumber
		
		
///	table a5: impact on endline assessments by awc attendance

	/// load data 1
	
		use ${data}irtssanaldata, clear

	///	panels a and b
		
		///	run analysis
		
			eststo clear
				
			foreach s in $subjects {
				qui eststo est3a_`s': areg stdf_`s'3 treat std_`s'1 ///
				awc_std_`s'1 $cov [aw=w] if (hh_sample==1| ///
				hh_sample==2) & awc_att==1, absorb(stratum) ///
				cluster(digitcode)						
				
				qui su stdf_`s'3 [aw=w] if treat==0 & ///
				(hh_sample==1|hh_sample==2) & awc_att==1
				qui estadd scalar ymean=r(mean)
				
				qui eststo est3b_`s': areg stdf_`s'3 treat std_`s'1 ///
				awc_std_`s'1 $cov [aw=w] if (hh_sample==1| ///
				hh_sample==2) & awc_att==0, absorb(stratum) ///
				cluster(digitcode)						
				
				qui su stdf_`s'3 [aw=w] if treat==0 & ///
				(hh_sample==1|hh_sample==2) & awc_att==0
				qui estadd scalar ymean=r(mean)				
			}	
		
		///	print table

			qui esttab est3a_* using ${tables}hetitt_byatt.tex, replace f ///
			refcat(treat ///
			"\emph{A. Attends AWC}\\\underline{HH assessments} (N=1129)", ///
			nolabel) nolines posthead(\midrule) rename(std_lang1 ///
			"std_math1" std_exec1 "std_math1" std_comb1 "std_math1") ///
			booktabs label se r2 drop(_cons std_* awc_* $cov) ///
			width(\hsize) b(%9.3f) nonotes star(* 0.1 ** 0.05 *** 0.01) ///
			mlabels("\shortstack{Math}" "\shortstack{Language}" ///
			"\shortstack{Executive\\ function}" ///
			"\shortstack{Composite\\ score}") stats(N ymean, fmt(%9.0f ///
			%9.3f) labels("N (children)" "Control mean"))
			
			qui esttab est3b_* using ${tables}hetitt_byatt.tex, append f ///
			refcat(treat ///
			"\emph{B. Does not attend}\\\underline{HH assessments} (N=946)", ///
			nolabel) nolines posthead(\midrule) rename(std_lang1 ///
			"std_math1" std_exec1 "std_math1" std_comb1 "std_math1") ///
			booktabs label se r2 drop(_cons std_* awc_* $cov) ///
			width(\hsize) b(%9.3f) nonotes star(* 0.1 ** 0.05 *** 0.01) ///
			stats(N ymean, fmt(%9.0f %9.3f) labels("N (children)" ///
			"Control mean")) mlabels(none) nonumber

/// table a.6: impact on endline WAZ/HAZ scores - common sample

	///	store weight globals
		
		global weight_awc	""
		global weight_hh	"[aw=w]"
	
	///	store title globals
	
		global title_waz	"WAZ score"
		global title_waz2	"Underweight\\(WAZ$<$-2)"
		global title_waz3	"Severely\\underweight\\(WAZ$<$-3)"
		global title_haz	"HAZ score"
		global title_haz2	"Stunted\\(HAZ$<$-2)"
		global title_haz3	"Severely\\stunted\\(HAZ$<$-3)"
		global title_whz	"WHZ score"
		global title_whz2	"Wasted\\(WHZ$<$-2)"
		global title_whz3	"Severely\\wasted\\(WHZ$<$-3)"

	/// load data 1
	
		use ${data}irtssanaldata, clear
	
	/// panel a
	
		///	run analysis 
		
			eststo clear	
			foreach v in waz06 low_waz very_low_waz {
				qui eststo esta_awc_`v': areg awc_`v' treat bas_waz06 ///
				$cov if bas_waz06!=. & awc_waz06!=. & hh_waz06!=. & ///
				(hh_sample==1|hh_sample==2), ///
				absorb(stratum) cluster(digitcode) 
					
				qui su awc_`v' if treat==0 & bas_waz06!=. & awc_waz06!=. ///
				& hh_waz06!=. & (hh_sample==1|hh_sample==2)
				qui estadd scalar ymean=r(mean)
					
				qui eststo esta_hh_`v': areg hh_`v' treat bas_waz06 ///
				$cov if bas_waz06!=. & awc_waz06!=. & hh_waz06!=. & ///
				(hh_sample==1|hh_sample==2), ///
				absorb(stratum) cluster(digitcode)
				
				qui su hh_`v' if treat==0 & bas_waz06!=. & awc_waz06!=. ///
				& hh_waz06!=. & (hh_sample==1|hh_sample==2)
				qui estadd scalar ymean=r(mean)
			}	
			
		///	print table
		
			qui esttab esta_awc_* using ${tables}estitt_nut_common.tex, ///
			replace f refcat(treat "\underline{AWC measurements} (N=790)", ///
			nolabel) nolines posthead(\midrule) booktabs label se ///
			drop(_cons $cov bas*) width(\hsize) b(%9.3f) nonotes ///
			star(* 0.1 ** 0.05 *** 0.01) stats(N ymean, fmt(%9.0f %9.3f ///
			%9.3f) labels("N (children)" "Control mean")) ///
			mlabels("\shortstack{${title_waz}}" ///
			"\shortstack{${title_waz2}}" "\shortstack{${title_waz3}}") 

			qui esttab esta_hh_* using ${tables}estitt_nut_common.tex, ///
			append f refcat(treat "\underline{HH measurements} (N=790)", ///
			nolabel) nolines booktabs label se drop(_cons $cov bas*) ///
			width(\hsize) b(%9.3f) nonotes star(* 0.1 ** 0.05 *** 0.01) ///
			stats(N ymean, fmt(%9.0f %9.3f %9.3f) labels("N (children)" ///
			"Control mean")) mlabels(none) nonumbers

	///	p-values
		
		/// prep data

			g in_bl_waz=bas_waz06!=. 
			g in_elawc_waz=awc_waz06!=. 
			g in_elhh_waz=hh_waz06!=. 
			
			gen i=_n
			expand 2
			bys i: gen s1=_n==1

			g weight=w if s1==1
			replace weight=w if s1==0
			
			foreach v in waz06 {
				g end_`v'=awc_`v' if s1==1
				replace end_`v'=hh_`v' if s1==0
					
				g bas2_`v'=bas_`v'*s1
				g bas3_`v'=bas_`v'*(1-s1)
			}
			
			foreach v in low_waz very_low_waz {
				g end_`v'=awc_`v' if s1==1
				replace end_`v'=hh_`v' if s1==0
			}

			g aww_pass102=aww_pass10*s1
			g aww_pass103=aww_pass10*(1-s1)

			g m_pass102=m_pass10*s1
			g m_pass103=m_pass10*(1-s1)
			
			g aww_exper2=aww_exper*s1
			g aww_exper3=aww_exper*(1-s1)

			g m_exper2=m_exper*s1
			g m_exper3=m_exper*(1-s1)
			
			g treat2=treat*s1
			g treat3=treat*(1-s1)
				
			egen sample_stratum=group(s1 stratum)
	
		///	run analysis

			foreach v in waz06 low_waz very_low_waz {
				qui eststo p_`v': areg end_`v' s1 bas2_waz06 bas3_waz06 ///
				treat2 treat3 aww_exper2 aww_exper3 aww_pass102 ///
				aww_pass103 m_exper2 m_exper3 m_pass102 m_pass103 if ///
				in_bl_waz==1 & in_elawc_waz==1 & in_elhh_waz==1 & ///
				(hh_sample==1|hh_sample==2), ///
				absorb(sample_stratum) cluster(digitcode)
				
				qui test treat2=treat3
				qui estadd scalar pval=r(p)
			}
			
		///	print table
		
			qui esttab p_waz06 p_low_waz p_very_low_waz ///
			using ${tables}estitt_nut_common.tex, append f nolines ///
			posthead(\midrule) drop(_cons s1 bas* treat* aww* m_*) ///
			booktabs label se width(\hsize) b(%9.3f) nonotes star(* 0.1 ** ///
			0.05 *** 0.01) stats(pval, fmt(%9.3f) ///
			labels("P-value (AWC $=$ HH)")) mlabels(none) nonumber

	/// load data 1
	
		use ${data}irtssanaldata, clear
			
	///	panel b
			
		///	run analysis 
		
			eststo clear	
			foreach v in haz06 low_haz very_low_haz {
				qui eststo estb_awc_`v': areg awc_`v' treat bas_haz06 ///
				$cov if bas_haz06!=. & awc_haz06!=. & hh_haz06!=. & ///
				(hh_sample==1|hh_sample==2), ///
				absorb(stratum) cluster(digitcode) 
					
				qui su awc_`v' if treat==0 & bas_haz06!=. & awc_haz06!=. ///
				& hh_haz06!=. & (hh_sample==1|hh_sample==2)
				qui estadd scalar ymean=r(mean)
					
				qui eststo estb_hh_`v': areg hh_`v' treat bas_haz06 ///
				$cov if bas_haz06!=. & awc_haz06!=. & hh_haz06!=. & ///
				(hh_sample==1|hh_sample==2), ///
				absorb(stratum) cluster(digitcode)
				
				qui su hh_`v' if treat==0 & bas_haz06!=. & awc_haz06!=. ///
				& hh_haz06!=. & (hh_sample==1|hh_sample==2)
				qui estadd scalar ymean=r(mean)
			}	
		
		///	print table
		
			qui esttab estb_awc_* using ${tables}estitt_nut_common.tex, ///
			append f refcat(treat "\underline{AWC measurements} (N=724)", ///
			nolabel) nolines posthead(\midrule) booktabs label se ///
			drop(_cons $cov bas*) width(\hsize) b(%9.3f) nonotes ///
			star(* 0.1 ** 0.05 *** 0.01) stats(N ymean, fmt(%9.0f %9.3f ///
			%9.3f) labels("N (children)" "Control mean")) ///
			mlabels("\shortstack{${title_haz}}" ///
			"\shortstack{${title_haz2}}" "\shortstack{${title_haz3}}") ///
			prehead(\toprule)

			qui esttab estb_hh_* using ${tables}estitt_nut_common.tex, ///
			append f refcat(treat "\underline{HH measurements} (N=724)", ///
			nolabel) nolines booktabs label se drop(_cons $cov bas*) ///
			width(\hsize) b(%9.3f) nonotes star(* 0.1 ** 0.05 *** 0.01) ///
			stats(N ymean, fmt(%9.0f %9.3f %9.3f) labels("N (children)" ///
			"Control mean")) mlabels(none) 

	///	p-values
		
		/// prep data

			g in_bl_haz=bas_haz06!=. 
			g in_elawc_haz=awc_haz06!=. 
			g in_elhh_haz=hh_haz06!=. 
			
			gen i=_n
			expand 2
			bys i: gen s1=_n==1

			g weight=w if s1==1
			replace weight=w if s1==0
			
			foreach v in haz06 {
				g end_`v'=awc_`v' if s1==1
				replace end_`v'=hh_`v' if s1==0
					
				g bas2_`v'=bas_`v'*s1
				g bas3_`v'=bas_`v'*(1-s1)
			}
			
			foreach v in low_haz very_low_haz {
				g end_`v'=awc_`v' if s1==1
				replace end_`v'=hh_`v' if s1==0
			}

			g aww_pass102=aww_pass10*s1
			g aww_pass103=aww_pass10*(1-s1)

			g m_pass102=m_pass10*s1
			g m_pass103=m_pass10*(1-s1)
			
			g aww_exper2=aww_exper*s1
			g aww_exper3=aww_exper*(1-s1)

			g m_exper2=m_exper*s1
			g m_exper3=m_exper*(1-s1)
			
			g treat2=treat*s1
			g treat3=treat*(1-s1)
				
			egen sample_stratum=group(s1 stratum)
	
		///	run analysis

			foreach v in haz06 low_haz very_low_haz {
				qui eststo p_`v': areg end_`v' s1 bas2_haz06 bas3_haz06 ///
				treat2 treat3 aww_exper2 aww_exper3 aww_pass102 ///
				aww_pass103 m_exper2 m_exper3 m_pass102 m_pass103 if ///
				in_bl_haz==1 & in_elawc_haz==1 & in_elhh_haz==1 & ///
				(hh_sample==1|hh_sample==2), ///
				absorb(sample_stratum) cluster(digitcode)
				
				qui test treat2=treat3
				qui estadd scalar pval=r(p)
			}
			
		///	print table
		
			qui esttab p_haz06 p_low_haz p_very_low_haz ///
			using ${tables}estitt_nut_common.tex, append f nolines ///
			posthead(\midrule) drop(_cons s1 bas* treat* aww* m_*) ///
			booktabs label se width(\hsize) b(%9.3f) nonotes star(* 0.1 ** ///
			0.05 *** 0.01) stats(pval, fmt(%9.3f) ///
			labels("P-value (AWC $=$ HH)")) mlabels(none) nonumber	

/// table a7: ipw-weighted impact on endline WAZ/HAZ scores
		
	///	store weight globals
		
		global weight_awc	""
		global weight_hh	"[aw=w]"
	
	///	store title globals
	
		global title_waz	"WAZ score"
		global title_waz2	"Underweight\\(WAZ$<$-2)"
		global title_waz3	"Severely\\underweight\\(WAZ$<$-3)"
		global title_haz	"HAZ score"
		global title_haz2	"Stunted\\(HAZ$<$-2)"
		global title_haz3	"Severely\\stunted\\(HAZ$<$-3)"
		global title_whz	"WHZ score"
		global title_whz2	"Wasted\\(WHZ$<$-2)"
		global title_whz3	"Severely\\wasted\\(WHZ$<$-3)"
				
	/// load data 1
	
		use ${data}irtssanaldata, clear
	
	/// est prob of followup
	
		///	mean-inpute for children w/o ach w/waz or haz
		
			foreach s in $subjects {
				replace std_`s'1=0 if  std_`s'1==. & bas_waz06!=.
				*ajg: this is to ensure ns are comparable
			}
	
		/// awc followup

			gen inv_prob2=.
			
			foreach t in 0 1 {
				probit fol_elawc female bas_waz06 bas_haz06 ///
				i.stratum if treat==`t', robust
				
				predict prob2`t' if treat==`t'				
				
				replace inv_prob2=1/prob2`t' if treat==`t'
			}
		
		///	hh followup
	
			gen inv_prob3=.
				
			foreach t in 0 1 {	
				probit fol_elhh female bas_waz06 bas_haz06 ///
				i.stratum if treat==`t', robust
			
				predict prob3`t' if treat==`t'
			
				replace inv_prob3=1/prob3`t' if treat==`t'
			}
			
			g combwt=w*inv_prob3
	
	/// panel a
	
		///	run analysis 
		
			eststo clear	
			foreach v in waz06 low_waz very_low_waz {
				qui eststo awc_`v': areg awc_`v' treat bas_waz06 $cov ///
				[pw=inv_prob2], absorb(stratum) cluster(digitcode)
						
				qui su awc_`v' if treat==0 & bas_waz06!=.
				qui estadd scalar ymean=r(mean)
				
				qui eststo hh_`v': areg hh_`v' treat bas_waz06 $cov ///
				[pw=combwt] if (hh_sample==1| ///
				hh_sample==2), absorb(stratum) cluster(digitcode)
						
				qui su hh_`v' [aw=w] if treat==0 & bas_waz06!=. & ///
				(hh_sample==1|hh_sample==2)
				qui estadd scalar ymean=r(mean)	
			}	
			
		///	print table
		
			qui esttab awc_* using ${tables}estitt_nut_ipw.tex, replace f ///
			refcat(treat "\underline{AWC measurements} (N=1513)", ///
			nolabel) nolines posthead(\midrule) booktabs label se ///
			drop(_cons $cov bas*) width(\hsize) b(%9.3f) nonotes ///
			star(* 0.1 ** 0.05 *** 0.01) stats(ymean, fmt(%9.3f) ///
			labels("Control mean")) mlabels("\shortstack{${title_waz}}" ///
			"\shortstack{${title_waz2}}" "\shortstack{${title_waz3}}") 

			qui esttab hh_* using ${tables}estitt_nut_ipw.tex, append f ///
			refcat(treat "\underline{HH measurements} (N=1997)", ///
			nolabel) nolines booktabs label se drop(_cons $cov bas*) ///
			width(\hsize) b(%9.3f) nonotes star(* 0.1 ** 0.05 *** 0.01) ///
			stats(ymean, fmt(%9.3f) labels("Control mean")) mlabels(none) ///
			nonumbers

	/// load data 1
	
		use ${data}irtssanaldata, clear
	
	/// est prob of followup
	
		///	mean-inpute for children w/o ach w/waz or haz
		
			foreach s in $subjects {
				replace std_`s'1=0 if  std_`s'1==. & bas_haz06!=.
				*ajg: this is to ensure ns are comparable
			}
	
		/// awc followup

			gen inv_prob2=.
			
			foreach t in 0 1 {
				probit fol_elawc female bas_waz06 bas_haz06 ///
				i.stratum if treat==`t', robust
				
				predict prob2`t' if treat==`t'				
				
				replace inv_prob2=1/prob2`t' if treat==`t'
			}
		
		///	hh followup
	
			gen inv_prob3=.
				
			foreach t in 0 1 {	
				probit fol_elhh female bas_waz06 bas_haz06 ///
				i.stratum if treat==`t', robust
			
				predict prob3`t' if treat==`t'
			
				replace inv_prob3=1/prob3`t' if treat==`t'
			}
			
			g combwt=w*inv_prob3
	
	/// panel a
	
		///	run analysis 
		
			eststo clear	
			foreach v in haz06 low_haz very_low_haz {
				qui eststo awc_`v': areg awc_`v' treat bas_haz06 $cov ///
				[pw=inv_prob2], absorb(stratum) cluster(digitcode)
						
				qui su awc_`v' if treat==0 & bas_haz06!=.
				qui estadd scalar ymean=r(mean)
				
				qui eststo hh_`v': areg hh_`v' treat bas_haz06 $cov ///
				[pw=combwt] if (hh_sample==1| ///
				hh_sample==2), absorb(stratum) cluster(digitcode)
						
				qui su hh_`v' [aw=w] if treat==0 & bas_haz06!=. & ///
				(hh_sample==1|hh_sample==2)
				qui estadd scalar ymean=r(mean)	
			}	
			
		///	print table
		
			qui esttab awc_* using ${tables}estitt_nut_ipw.tex, append f ///
			refcat(treat ///
			"\underline{AWC measurements} (N=1388)", ///
			nolabel) nolines posthead(\midrule) booktabs label se ///
			drop(_cons $cov bas*) width(\hsize) b(%9.3f) nonotes ///
			star(* 0.1 ** 0.05 *** 0.01) stats(ymean, fmt(%9.3f) ///
			labels("Control mean")) mlabels("\shortstack{${title_haz}}" ///
			"\shortstack{${title_haz2}}" "\shortstack{${title_haz3}}") ///
			prehead(\toprule)

			qui esttab hh_* using ${tables}estitt_nut_ipw.tex, append f ///
			refcat(treat "\underline{HH measurements} (N=1988)", ///
			nolabel) nolines booktabs label se drop(_cons $cov bas*) ///
			width(\hsize) b(%9.3f) nonotes star(* 0.1 ** 0.05 *** 0.01) ///
			stats(ymean, fmt(%9.3f) labels("Control mean")) mlabels(none) ///
			nonumbers
			
/// tables a8-a9: heterogeneous impact on endline WAZ/HAZ scores
	
	/// load data 1
	
		use ${data}irtssanaldata, clear

	///	lab coeffs
	
		lab var treat		"Treatment"		
		lab var waz			"Covariate"
		lab var	haz			"Covariate"
		lab var i_waz		"Interaction"
		lab var i_haz		"Interaction"
		*ajg: i do this to place all bl coeffs on the same row
		
	foreach m in waz haz {
		
	/// panel a
	
		///	run analysis 
		
			eststo clear
			
			qui eststo awc_`m': areg awc_`m'06 treat `m' i_`m' ///
			$cov, absorb(stratum) cluster(digitcode)
			
			qui eststo hh_`m': areg hh_`m'06 treat `m' i_`m' ///
			$cov [aw=w] if (hh_sample==1|hh_sample==2), ///
			absorb(stratum) cluster(digitcode)
			
			foreach i in female months moth10 vacancy children comb {
				qui eststo awc_`m'_`i': areg awc_`m'06 treat bas_`m'06 ///
				`i' i_`i' $cov, absorb(stratum) cluster(digitcode)				
				
				qui eststo hh_`m'_`i': areg hh_`m'06 treat bas_`m'06 ///
				`i' i_`i' $cov [aw=w] if (hh_sample==1| ///
				hh_sample==2), absorb(stratum) cluster(digitcode)
			}	
			/*ajg: here, i include comb in the loop because it is always
			the main effect; it is never a control variable */
			
		///	print table
		
			qui esttab awc_* using ${tables}hetitt_`m'.tex, replace f ///
			nolines posthead("\midrule" "\multicolumn{2}{l}{" ///
			"\emph{A. Complete sample}}\\" "\multicolumn{2}{l}{" ///
			"\underline{AWC measurements}}\\\addlinespace") ///
			booktabs label se b(%9.3f) nonotes star(* 0.1 ** 0.05 *** ///
			0.01) stats(N, fmt(%9.0f) labels("N (children)")) ///
			width(\hsize) drop(_cons $cov bas_*) rename(female "`m'" ///
			months "`m'" moth10 "`m'" vacancy "`m'" children "`m'" comb ///
			"`m'" i_female "i_`m'" i_months "i_`m'" i_moth10 "i_`m'" ///
			i_vacancy "i_`m'" i_children "i_`m'" i_comb "i_`m'") ///
			mlabels("\shortstack{Baseline\\ ${`m'lab}}" "Female" ///
			"\shortstack{Age\\(months)}" "\shortstack{Mother\\education}" ///
			"\shortstack{AWW\\vacancy}" "\shortstack{Class\\size}" ///
			"\shortstack{Baseline\\score}")
			
			qui esttab hh_* using ${tables}hetitt_`m'.tex, append f ///
			nolines posthead("\addlinespace \multicolumn{2}{l}{" ///
			"\underline{HH measurements}}\\\addlinespace") ///
			booktabs label se b(%9.3f) nonotes star(* 0.1 ** 0.05 *** ///
			0.01) stats(N, fmt(%9.0f) labels("N (children)")) ///
			width(\hsize) drop(_cons $cov bas_*) rename(female "`m'" ///
			months "`m'" moth10 "`m'" vacancy "`m'" children "`m'" comb ///
			"`m'" i_female "i_`m'" i_months "i_`m'" i_moth10 "i_`m'" ///
			i_vacancy "i_`m'" i_children "i_`m'" i_comb "i_`m'") ///
			mlabels(none) nonumber
			
	///	panel b		
	
		/// run analysis
		
			eststo clear
			
			qui eststo awc_`m': areg awc_`m'06 treat `m' i_`m' ///
			$cov if bas_`m'06!=. & awc_`m'06!=. & hh_`m'06!=. & ///
			(hh_sample==1|hh_sample==2), ///
			absorb(stratum) cluster(digitcode)
			
			qui eststo hh_`m': areg hh_`m'06 treat `m' i_`m' ///
			$cov [aw=w] if bas_`m'06!=. & awc_`m'06!=. & hh_`m'06!=. & ///
			(hh_sample==1|hh_sample==2), ///
			absorb(stratum) cluster(digitcode)	
			
			foreach i in female months moth10 vacancy children comb {
				qui eststo awc_`m'_`i': areg awc_`m'06 treat bas_`m'06 ///
				`i' i_`i' $cov if bas_`m'06!=. & awc_`m'06!=. & ///
				hh_`m'06!=. & (hh_sample==1|hh_sample==2), ///
				absorb(stratum) cluster(digitcode)				
				
				qui eststo hh_`m'_`i': areg hh_`m'06 treat bas_`m'06 ///
				`i' i_`i' $cov if bas_`m'06!=. & awc_`m'06!=. & ///
				hh_`m'06!=. & (hh_sample==1|hh_sample==2), ///
				absorb(stratum) cluster(digitcode)
			}
		
		///	print table
	
			qui esttab awc_* using ${tables}hetitt_`m'.tex, append f ///
			nolines posthead("\midrule" "\multicolumn{2}{l}{" ///
			"\emph{B. Common sample}}\\" "\multicolumn{2}{l}{" ///
			"\underline{AWC measurements}}\\\addlinespace") booktabs label ///
			se width(\hsize) b(%9.3f) nonotes star(* 0.1 ** 0.05 *** 0.01) ///
			drop(_cons $cov bas_*) rename(female "`m'" months ///
			"`m'" moth10 "`m'" vacancy "`m'" children "`m'" comb "`m'" ///
			i_female "i_`m'" i_months "i_`m'" i_moth10 "i_`m'" i_vacancy ///
			"i_`m'" i_children "i_`m'" i_comb "i_`m'") ///
			mlabels(none) nonumbers noobs

			qui esttab hh_* using ${tables}hetitt_`m'.tex, append f ///
			nolines posthead("\addlinespace \multicolumn{2}{l}{" ///
			"\underline{HH measurements}}\\\addlinespace") booktabs label ///
			se b(%9.3f) nonotes star(* 0.1 ** 0.05 *** 0.01) ///
			stats(N, fmt(%9.0f) labels("N (children)")) width(\hsize) ///
			drop(_cons $cov bas_*) rename(female "`m'" months "`m'" moth10 ///
			"`m'" vacancy "`m'" children "`m'" comb "`m'" i_female "i_`m'" ///
			i_months "i_`m'" i_moth10 "i_`m'" i_vacancy "i_`m'" i_children ///
			"i_`m'" i_comb "i_`m'") mlabels(none) nonumber
	}

///	tables a10, a11: impact on endline HAZ/WAZ scores by awc attendance

	/// load data 1
	
		use ${data}irtssanaldata, clear

	foreach m in waz haz {
	
	///	panel a
	
		/// run analysis
		
			eststo clear
	
			foreach v in `m'06 low_`m' very_low_`m' {
				qui eststo esta_hh_`v': areg hh_`v' treat bas_`m'06 $cov ///
				[aw=w] if awc_att==1 & (hh_sample==1| ///
				hh_sample==2), absorb(stratum) cluster(digitcode)
							
				qui su hh_`v' [aw=w] if treat==0 & bas_`m'06!=. & ///
				awc_att==1 & (hh_sample==1|hh_sample==2)
				qui estadd scalar ymean=r(mean)
			}	
			
		///	print table
		
			qui esttab esta_hh_* using ${tables}hetitt_`m'_byatt.tex, ///
			replace f refcat(treat ///
			"\emph{A. Attends AWC}\\\underline{HH measurements} (N=1081)", ///
			nolabel) nolines posthead(\midrule) booktabs label se ///
			drop(_cons $cov bas*) width(\hsize) b(%9.3f) nonotes ///
			star(* 0.1 ** 0.05 *** 0.01) stats(N ymean, fmt(%9.0f %9.3f ///
			%9.3f) labels("N (children)" "Control mean")) ///
			mlabels("\shortstack{${title_`m'}}" ///
			"\shortstack{${title_`m'2}}" "\shortstack{${title_`m'3}}")
			
	///	panel b
	
		/// run analysis
		
			foreach v in `m'06 low_`m' very_low_`m' {				
				qui eststo estb_hh_`v': areg hh_`v' treat bas_`m'06 ///
				$cov [aw=w] if awc_att==0 & (hh_sample==1| ///
				hh_sample==2), absorb(stratum) cluster(digitcode)
						
				qui su hh_`v' [aw=w] if treat==0 & bas_`m'06!=. & ///
				awc_att==0 & (hh_sample==1|hh_sample==2)
				qui estadd scalar ymean=r(mean)
			}
			
		///	print table
		
			qui esttab estb_hh_* using ${tables}hetitt_`m'_byatt.tex, ///
			append f refcat(treat ///
			"\emph{B. Does not attend}\\\underline{HH measurements} (N=909)", ///
			nolabel) nolines posthead(\midrule) booktabs label se ///
			drop(_cons $cov bas*) width(\hsize) b(%9.3f) nonotes ///
			star(* 0.1 ** 0.05 *** 0.01) stats(N ymean, fmt(%9.0f %9.3f ///
			%9.3f) labels("N (children)" "Control mean")) mlabels(none) ///
			nonumbers		
	}	

///	table a12: regression-based association between nutrition status and learning outcomes	

	/// load data 1
	
		use ${data}irtssanaldata, clear
	
	///	panel a
	
		/// run analysis
		
			foreach s in $subjects {
				foreach m in waz haz {
					foreach v in `m'06 low_`m' very_low_`m' {
						reg std_`s'1 bas_`v'
					}	
				}	
			}	
	
///	table b1: comparison of in- and out-of-sample districts


///	variable globals
	
	# delimit ;
 
	global	dist_vars		"prop_normal prop_muw prop_suw prop_ow
							districtpopulation";
							
	# delimit cr	

	global	distpop 		${data}tn_districts_wpop_20160415.xlsx
	global	distnut			${data}tn_districts_wnut_20160415.xlsx

/// import population data
	*ajg: sampled districts are 604, 613, 614, 625

	import excel $distpop, sheet("Sheet1") clear firstrow

	/// make vars lower case
	
		rename *, lower
	
	/// drop chennai
	
		drop if district=="Chennai"
		
	/// en district names
	
		en district, g(district_name)
		/*ajg: i can't figure out why i can't match the var "district" from
		this dataset and the var "districtname" from the other dataset;
		i thought it was an extra space or format but neither worked */
		drop district
		
	///	save as temp file
		
		tempfile temp_dist_pop
		save `temp_dist_pop'	
		
///	import nutrition data
	
	import excel $distnut, sheet("Sheet1") clear firstrow	
	
	/// make vars lower case
	
		rename *, lower
	
	/// drop chennai
	
		drop if districtname=="Chennai"
	
	/// en district names
	
		en districtname, g(district_name)
		/*ajg: i can't figure out why i can't match the var "district" from
		this dataset and the var "districtname" from the other dataset;
		i thought it was an extra space or format but neither worked */
		drop districtname
	
	///	merge w/pop data
	
		mer 1:1 district_name using `temp_dist_pop', nogen
		/*-----------------------------------------
		Not matched                             0
		Matched                                31  
		-----------------------------------------*/

	/// g dummy for in-sample districts
	
		g sample=(did==604|did==613|did==614|did==625)
		*ajg: sampled districts are 604, 613, 614, 625

	///	g prop vars
		*ajg: can't compare #s; need to compare proportions
		
		foreach v in normal muw suw ow {
			g prop_`v'=`v'/total
			drop `v'
		}
		
	///	recalculate pop figures
	
		replace districtpopulation=districtpopulation/1000
	
	///	label vars
	
		lab var	prop_normal	"Proportion of children with normal weight"
		lab var	prop_muw	"Proportion of moderately underweight children"
		lab var	prop_suw	"Proportion of severely underweight chidlren"
		lab var	prop_ow		"Proportion of overweight children"
		lab var	districtpop	"Number of residents in the district (in 1000s)"
	
	///	table b13: comparison of in- and out-of-sample districts	
	
		balancetable sample $dist_vars using ${tables}checkextval.tex, ///
		replace booktabs leftctitle(" ") ///
		ctit("Out-of-sample" "In-sample" "Difference" "N") ///
		observationscolumn noobs varlabels sd(par([ ]))

///	table b2: facilitator self-reports on time allocation from intervention monitoring

	/// load data 1
	
		use ${data}imanaldata1, clear
	
	///	print table
	
		balancetable (mean) $intmon_vars2a ///
		using ${tables}sustats_intmon2.tex, replace booktabs varlabels ///
		noobs nolines nofoot sd(par([ ])) ///
		prehead("\begin{tabular}{l*{1}c}\toprule Hours spent last week...") ///
		posthead("\midrule \emph{A. Education}\\") 

		balancetable (mean)	$intmon_vars2b ///
		using ${tables}sustats_intmon2.tex, append booktabs varlabels noobs ///
		nolines nohead nofoot sd(par([ ])) prehead("\midrule") ///
		posthead("\emph{B. Health and nutrition}\\") 

		balancetable (mean) $intmon_vars2c ///
		using ${tables}sustats_intmon2.tex, append booktabs varlabels noobs ///
		nolines nohead nofoot sd(par([ ])) prehead("\midrule") ///
		posthead("\emph{C. Administrative}\\") 

		balancetable (mean) $intmon_vars2d ///
		using ${tables}sustats_intmon2.tex, append booktabs varlabels ///
		nolines nohead sd(par([ ])) postfoot("\bottomrule\end{tabular}") ///
		prehead("\midrule") leftcobservations("N (centers)")		
					
///	table b3: impact on endline assessments (proportion-correct scores)

	/// load data 1
	
		use ${data}irtssanaldata, clear
	
	/// panel a
	
		/// run analysis

			eststo clear
			
			foreach s in $subjects {
				eststo est2_`s': areg per_`s'2 treat per_`s'1 ///
				awc_per_`s'1 $cov, absorb(stratum) cluster(digitcode)
				
				qui su per_`s'2 if treat==0 & per_`s'1!=.
				qui estadd scalar ymean=r(mean)
				
				eststo est3_`s': areg per_`s'3 treat per_`s'1 ///
				awc_per_`s'1 $cov [aw=w] if (hh_sample==1| ///
				hh_sample==2), absorb(stratum) cluster(digitcode)
				
				qui su per_`s'3 [aw=w] if treat==0 & per_`s'1!=. & ///
				(hh_sample==1|hh_sample==2)
				qui estadd scalar ymean=r(mean)				
			}
			
		///	print table
		
			qui esttab est2_* using ${tables}estitt_per.tex, replace f ///
			refcat(treat ///
			"\emph{A. Complete sample}\\\underline{AWC assessments} (N=1514)", ///
			nolabel) nolines posthead(\midrule) rename(per_lang1 ///
			"per_math1" per_exec1 "per_math1" per_comb1 "per_math1") ///
			booktabs label se r2 drop(_cons per_* awc_* $cov) ///
			width(\hsize) b(%9.3f) nonotes star(* 0.1 ** 0.05 *** 0.01) ///
			mlabels("\shortstack{Math}" "\shortstack{Language}" ///
			"\shortstack{Executive\\ function}" ///
			"\shortstack{Composite\\ score}") stats(ymean, fmt(%9.3f) ///
			labels("Control mean"))
		
			qui esttab est3_* using ${tables}estitt_per.tex, append f ///
			refcat(treat "\underline{HH assessments} (N=2075)", nolabel) ///
			nolines rename(per_lang1 "per_math1" per_exec1 "per_math1" ///
			per_comb1 "per_math1") booktabs label se r2 drop(_cons ///
			per_* awc_* $cov) width(\hsize) b(%9.3f) nonotes ///
			star(* 0.1 ** 0.05 *** 0.01) stats(ymean, fmt(%9.3f) ///
			labels("Control mean")) mlabels(none) nonumber
	
	///	panel b
	
		/// run analysis
	
			eststo clear
			
			foreach s in $subjects {
				qui eststo est2_`s': areg per_`s'2 treat per_`s'1 ///
				awc_per_`s'1 $cov if in_elawc==1 & in_elhh==1 & ///
				(hh_sample==1|hh_sample==2), ///
				absorb(stratum) cluster(digitcode)
	
				qui su per_`s'2 if treat==0 & per_`s'1!=. & in_elawc==1 & ///
				in_elhh==1 & (hh_sample==1|hh_sample==2)
				qui estadd scalar ymean=r(mean)
	
				qui eststo est3_`s': areg per_`s'3 treat per_`s'1 ///
				awc_per_`s'1 $cov if in_elawc==1 & in_elhh==1 & ///
				(hh_sample==1|hh_sample==2), ///
				absorb(stratum) cluster(digitcode)				

				qui su per_`s'3 if treat==0 & per_`s'1!=. & in_elawc==1 & ///
				in_elhh==1 & (hh_sample==1|hh_sample==2)
				qui estadd scalar ymean=r(mean)				
			}
			
		///	print table
		
			qui esttab est2_* using ${tables}estitt_per.tex, append f ///
			refcat(treat ///
			"\emph{B. Common sample}\\\underline{AWC assessments} (N=791)", ///
			nolabel) nolines posthead(\midrule) rename(per_lang1 ///
			"per_math1" per_exec1 "per_math1" per_comb1 "per_math1") ///
			booktabs label se drop(_cons per_* awc_* $cov) width(\hsize) ///
			b(%9.3f) nonotes star(* 0.1 ** 0.05 *** 0.01) stats(ymean, ///
			fmt(%9.3f) labels("Control mean")) mlabels(none) nonumber

			qui esttab est3_* using ${tables}estitt_per.tex, append f ///
			refcat(treat "\underline{HH assessments} (N=791)", nolabel) ///
			nolines rename(per_lang1 "per_math1" per_exec1 "per_math1" ///
			per_comb1 "per_math1") booktabs label se drop(_cons per_* ///
			awc_* $cov) width(\hsize) b(%9.3f) nonotes star(* 0.1 ** ///
			0.05 *** 0.01) stats(ymean, fmt(%9.3f) labels("Control mean")) ///
			mlabels(none) nonumber
	
	///	p-values
		
		/// prep data
				
			gen i=_n
			expand 2
			bys i: gen s1=_n==1
				
			foreach s in $subjects {
				g per_`s'=per_`s'2 if s1==1
				replace per_`s'=per_`s'3 if s1==0
						
				g bas2_per_`s'=per_`s'1*s1
				g bas3_per_`s'=per_`s'1*(1-s1)
					
				g awc2_per_`s'=awc_per_`s'1*s1
				g awc3_per_`s'=awc_per_`s'1*(1-s1)				
			}
						
			g weight=1 if s1==1
			replace weight=w if s1==0

			g aww_pass102=aww_pass10*s1
			g aww_pass103=aww_pass10*(1-s1)

			g m_pass102=m_pass10*s1
			g m_pass103=m_pass10*(1-s1)
				
			g aww_exper2=aww_exper*s1
			g aww_exper3=aww_exper*(1-s1)

			g m_exper2=m_exper*s1
			g m_exper3=m_exper*(1-s1)
				
			g treat2=treat*s1
			g treat3=treat*(1-s1)
					
			egen sample_stratum=group(s1 stratum)
		
		///	run analysis
				
			foreach s in $subjects {
				qui eststo pper_`s': areg per_`s' s1 bas2_per_`s' ///
				bas3_per_`s' awc2_per_`s' awc3_per_`s' treat2 treat3 ///
				aww_exper2 aww_exper3 m_exper2 m_exper3 aww_pass102 ///
				aww_pass103 m_pass102 m_pass103 [aw=weight] if ///
				in_elawc==1 & in_elhh==1 & (hh_sample==1| ///
				hh_sample==2), absorb(sample_stratum) ///
				cluster(digitcode)
						
				qui test treat2=treat3
				qui estadd scalar pval=r(p)
			}	
			
		///	print table

			qui esttab pper_* using ${tables}estitt_per.tex, append f ///
			nolines posthead(\midrule) drop(_cons s1 bas* awc* treat* ///
			aww* m_*) booktabs label se width(\hsize) b(%9.3f) nonotes ///
			star(* 0.1 ** 0.05 *** 0.01) stats(pval, fmt(%9.3f) ///
			labels("P-value (AWC $=$ HH)")) mlabels(none) nonumber	

/// table b4: impact on endline assessments (any-correct answers)

	/// load data 1
	
		use ${data}irtssanaldata, clear
		
	/// panel a
	
		/// run analysis

			eststo clear
			
			foreach s in $subjects {
				eststo est2_`s': areg any_`s'2 treat any_`s'1 ///
				awc_std_`s'1 $cov, absorb(stratum) cluster(digitcode)
				
				su any_`s'2 if treat==0 & any_`s'1!=.
				qui estadd scalar ymean=r(mean)
				
				eststo est3_`s': areg any_`s'3 treat any_`s'1 ///
				awc_std_`s'1 $cov [aw=w] if (hh_sample==1| ///
				hh_sample==2), absorb(stratum) cluster(digitcode)
				
				su any_`s'3 [aw=w] if treat==0 & any_`s'1!=. & ///
				(hh_sample==1|hh_sample==2)
				qui estadd scalar ymean=r(mean)
			}
			
		///	print table
		
			qui esttab est2_* using ${tables}estitt_any.tex, replace f ///
			refcat(treat ///
			"\emph{A. Complete sample}\\\underline{AWC assessments} (N=1514)", ///
			nolabel) nolines posthead(\midrule) rename(std_lang1 ///
			"std_math1" std_exec1 "std_math1" std_comb1 "std_math1") ///
			booktabs label se r2 drop(_cons any_* awc_* $cov) ///
			width(\hsize) b(%9.3f) nonotes star(* 0.1 ** 0.05 *** 0.01) ///
			mlabels("\shortstack{Math}" "\shortstack{Language}" ///
			"\shortstack{Executive\\ function}" ///
			"\shortstack{Composite\\ score}") ///
			stats(ymean, fmt(%9.3f) labels("Control mean"))
		
			qui esttab est3_* using ${tables}estitt_any.tex, append f ///
			refcat(treat "\underline{HH assessments} (N=2075)", nolabel) ///
			nolines rename(std_lang1 "std_math1" std_exec1 "std_math1" ///
			std_comb1 "std_math1") booktabs label se r2 drop(_cons ///
			any_* awc_* $cov) width(\hsize) b(%9.3f) nonotes ///
			star(* 0.1 ** 0.05 *** 0.01) stats(ymean, fmt(%9.3f) ///
			labels("Control mean")) mlabels(none) nonumber
	
	///	panel b
	
		/// run analysis
	
			eststo clear
			
			foreach s in $subjects {
				qui eststo est2_`s': areg any_`s'2 treat any_`s'1 ///
				awc_std_`s'1 $cov if in_elawc==1 & in_elhh==1 & ///
				(hh_sample==1|hh_sample==2), ///
				absorb(stratum) cluster(digitcode)
	
				qui su any_`s'2 if treat==0 & any_`s'1!=. & in_elawc==1 & ///
				in_elhh==1 & (hh_sample==1|hh_sample==2)
				qui estadd scalar ymean=r(mean)
	
				qui eststo est3_`s': areg any_`s'3 treat any_`s'1 ///
				awc_std_`s'1 $cov if in_elawc==1 & in_elhh==1 & ///
				(hh_sample==1|hh_sample==2), ///
				absorb(stratum) cluster(digitcode)				

				qui su any_`s'3 if treat==0 & any_`s'1!=. & in_elawc==1 & ///
				in_elhh==1 & (hh_sample==1|hh_sample==2)
				qui estadd scalar ymean=r(mean)				
			}
			
		///	print table
		
			qui esttab est2_* using ${tables}estitt_any.tex, append f ///
			refcat(treat ///
			"\emph{B. Common sample}\\\underline{AWC assessments} (N=791)", ///
			nolabel) nolines posthead(\midrule) rename(std_lang1 ///
			"std_math1" std_exec1 "std_math1" std_comb1 "std_math1") ///
			booktabs label se drop(_cons any_* awc_* $cov) width(\hsize) ///
			b(%9.3f) nonotes star(* 0.1 ** 0.05 *** 0.01) stats(ymean, ///
			fmt(%9.3f) labels("Control mean")) mlabels(none) nonumber

			qui esttab est3_* using ${tables}estitt_any.tex, append f ///
			refcat(treat "\underline{HH assessments} (N=791)", nolabel) ///
			nolines rename(std_lang1 "std_math1" std_exec1 "std_math1" ///
			std_comb1 "std_math1") booktabs label se drop(_cons any_* ///
			awc_* $cov) width(\hsize) b(%9.3f) nonotes star(* 0.1 ** ///
			0.05 *** 0.01) stats(ymean, fmt(%9.3f) labels("Control mean")) ///
			mlabels(none) nonumber	

///	table b5: impact on endline assessments on the same sample and items

	/// load data 1
	
		use ${data}irtssanaldata, clear
		
	///	panel a
	
		/// run analysis
	
			eststo clear
			
			foreach s in $subjects_nc {
				qui eststo `s'4_per: areg per_`s'4 treat per_`s'1 ///
				awc_per_`s'1 $cov if in_elawc==1 & in_elhh==1 & ///
				(hh_sample==1|hh_sample==2), ///
				absorb(stratum) cluster(digitcode)

				qui su per_`s'4 if treat==0 & per_`s'1!=. & in_elawc==1 & ///
				in_elhh==1 & (hh_sample==1|hh_sample==2) 
				qui estadd scalar ymean=r(mean)

				qui eststo `s'4_std: areg stdc_`s'4 treat std_`s'1 ///
				awc_std_`s'1 $cov if in_elawc==1 & in_elhh==1 & ///
				(hh_sample==1|hh_sample==2), ///
				absorb(stratum) cluster(digitcode)
				/*ajg: can't do a loop for per+std because, for std, i
				specifically want to use the stdc vars as outcomes, not std */

				qui su stdc_`s'4 if treat==0 & std_`s'1!=. & in_elawc==1 & ///
				in_elhh==1 & (hh_sample==1|hh_sample==2) 
				qui estadd scalar ymean=r(mean)
				
			}

			qui eststo comb4_std: areg stdc_comb4 treat std_comb1 ///
			awc_std_comb1 $cov if (in_elawc==1 & in_elhh==1 & ///
			(hh_sample==1|hh_sample==2)), absorb(stratum)  ///
			cluster(digitcode)
			/*ajg: can't include these in the loop because i want to exclude
			comb per from the table (it's meaningles) */
					
			qui su stdc_comb4 if treat==0 & std_comb1!=. & in_elawc==1 & ///
			in_elhh==1 & (hh_sample==1|hh_sample==2)
			qui estadd scalar ymean=r(mean)
			
		///	print table
		
			qui esttab _all using ${tables}estitt_sameboth.tex, replace f ///
			refcat(treat "\emph{A. AWC} (N=791)", nolabel) nolines ///
			posthead(\midrule) rename(per_lang1 "per_math1" per_exec1 ///
			"per_math1" per_comb1 "per_math1" std_math1 "per_math1" ///
			std_lang1 "per_math1" std_exec1 "per_math1" std_comb1 ///
			"per_math1") booktabs label se r2 keep(treat) width(\hsize) ///
			b(%9.3f) nonotes star(* 0.1 ** 0.05 *** 0.01) stats(ymean, ///
			fmt(%9.3f) labels("Control mean")) ///	
			mlabels("\shortstack{Raw}" "\shortstack{Std.}" ///
			"\shortstack{Raw}" "\shortstack{Std.}" ///
			"\shortstack{Raw}" "\shortstack{Std.}" "\shortstack{Std.}") ///
			mgroups("Math" "Language" "\shortstack{Executive\\ function}" ///
			"\shortstack{Composite\\ score}", ///
			pattern(1 0 1 0 1 0 1) prefix(\multicolumn{@span}{c}{) ///
			suffix(}) span erepeat(\cmidrule(lr){@span}))
		
	///	panel b
	
		/// run analysis
	
			eststo clear
			
			foreach s in $subjects_nc {
				qui eststo `s'5_per: areg per_`s'5 treat per_`s'1 ///
				awc_per_`s'1 $cov if in_elawc==1 & in_elhh==1 & ///
				(hh_sample==1|hh_sample==2), ///
				absorb(stratum) cluster(digitcode)
				
				qui su per_`s'5 if treat==0 & per_`s'1!=. & in_elawc==1 & ///
				in_elhh==1 & (hh_sample==1|hh_sample==2)
				qui estadd scalar ymean=r(mean)				
			
				qui eststo `s'5_std: areg stdc_`s'5 treat std_`s'1 ///
				awc_std_`s'1 $cov if in_elawc==1 & in_elhh==1 & ///
				(hh_sample==1|hh_sample==2), ///
				absorb(stratum) cluster(digitcode)
				
				qui su stdc_`s'5 if treat==0 & std_`s'1!=. & in_elawc==1 & ///
				in_elhh==1 & (hh_sample==1|hh_sample==2)
				qui estadd scalar ymean=r(mean)				
			
			}
	
			qui eststo comb5_std: areg std_comb5 treat std_comb1 ///
			awc_std_comb1 $cov if (in_elawc==1 & in_elhh==1 & ////
			(hh_sample==1|hh_sample==2)), absorb(stratum) cluster(digitcode)
					
			qui su std_comb5 if treat==0 & std_comb1!=. & in_elawc==1 & ///
			in_elhh==1 & (hh_sample==1|hh_sample==2)
			qui estadd scalar ymean=r(mean)
	
		///	print table
						
			qui esttab _all using ${tables}estitt_sameboth.tex, append f ///
			refcat(treat "\emph{B. HH} (N=791)", nolabel) nolines ///
			posthead(\midrule) rename(per_lang1 "per_math1" per_exec1 ///
			"per_math1" per_comb1 "per_math1" std_math1 "per_math1" ///
			std_lang1 "per_math1" std_exec1 "per_math1" std_comb1 ///
			"per_math1") booktabs label se r2 keep(treat) width(\hsize) ///
			b(%9.3f) nonotes star(* 0.1 ** 0.05 *** 0.01) stats(ymean, ///
			fmt(%9.3f) labels("Control mean")) ///	
			mlabels(none) nonumber

	///	p-values
		
		/// prep data
		
			foreach v in $subjects {
				forval r = 4/5 {
					replace std_`v'`r'=stdc_`v'`r'
					*ajg: i do this to make looping easier below
				}
			}

			keep per* std* in_elawc in_elhh treat w digitcode childid ///
			stratum $cov hh_sample
			
			foreach v in $subjects {
				foreach t in per std {
					g bas_`t'_`v'=`t'_`v'1
				}
			}
			
			reshape long per_math per_lang per_exec per_comb std_math ///
			std_lang std_exec std_comb, i(digitcode childid stratum w ///
			treat $cov hh_sample) j(round)
			
			g awc=round==2|round==4
			
			g i_awc=treat*awc	
	
		///	run analysis

			eststo clear
		
			foreach v in $subjects {
				foreach t in per std {
					qui eststo p`v'_`t': areg `t'_`v' treat awc i_awc ///
					bas_`t'_`v' $cov if ((round==4 | round==5) & ///
					in_elawc==1 & in_elhh==1 & (hh_sample==1| ///
					hh_sample==2)), absorb(stratum) cluster(digitcode)
					
					qui estadd scalar ///
					pval=(2*ttail(e(df_r),(abs(_b[i_awc]/_se[i_awc]))))
				}
			}
		
		///	print table
		
			qui esttab pmath_* plang_* pexec_* pcomb_std using ///
			${tables}estitt_sameboth.tex, append f nolines ///
			posthead(\midrule) drop(_cons awc bas_* i_awc aww* treat ///
			$cov) booktabs label se r2 width(\hsize) b(%9.3f) nonotes ///
			star(* 0.1 ** 0.05 *** 0.01) stats(pval, fmt(%9.3f) ///
			labels("P-value (T $\times$ AWC)")) mlabels(none) nonumber
		
///	table b6-b7: impacts on WAZ/HAZ scores after dropping child outliers

	/// load data 1
	
		use ${data}irtssanaldata, clear
	
	///	set table globals
	
		global	post4	"\midrule"
		global	post3	"\midrule"
		global	post2	"\midrule"
		global	post1	""
		*ajg: i do this to avoid double line at the end of table
	
	foreach m in waz haz {
	
		eststo clear
	
		foreach v in `m'06 low_`m' very_low_`m' {
			
		///	all panels
			
			/// run analysis
			
				forval r=1/5 {
					qui eststo p1_`v'_r`r': areg awc_`v' treat bas_`m'06 ///
					$cov [aw=w] if bas_`m'06!=.  & awc_`m'06!=. & ///
					cout_`m'_r`r'~=1, absorb(stratum) cluster(digitcode) 
			
					qui su awc_`v' if treat==0 & bas_`m'06!=. & ///
					awc_`m'06!=. & cout_`m'_r`r'~=1
					qui estadd scalar ymean=r(mean)	
				}
				
			///	print table
			
				qui esttab p1_*_r5 using ${tables}estitt_`m'_out.tex, ///
				replace f refcat(treat ///
				"\emph{Dropping residuals $>$ 0.5}", nolabel) nolines ///
				posthead(\midrule) booktabs label se drop(_cons $cov bas*) ///
				width(\hsize) b(%9.3f) nonotes star(* 0.1 ** 0.05 *** ///
				0.01) stats(N ymean, fmt(%9.0f %9.3f) ///
				labels("N (children)" "Control mean")) ///
				mlabels("\shortstack{${title_`m'}}" ///
				"\shortstack{${title_`m'2}}" "\shortstack{${title_`m'3}}") ///
				postfoot(\midrule)
				
				foreach r in 4 3 2 1 {
					qui esttab  p1_*_r`r' using ///
					${tables}estitt_`m'_out.tex, append f refcat(treat ///
					"\emph{Dropping residuals $>$ 0.`r'}", nolabel) ///
					nolines booktabs label se drop(_cons $cov bas*) ///
					width(\hsize) b(%9.3f) nonotes star(* 0.1 ** ///
					0.05 *** 0.01) stats(N ymean, fmt(%9.0f %9.3f %9.3f) ///
					labels("N (children)" "Control mean")) mlabels(none) ///
					nonumbers postfoot(${post`r'})		
				}
		}		
	}
	
///	table b8-b9: impacts on winsorized WAZ/HAZ scores

	/// load data 1
	
		use ${data}irtssanaldata, clear

	///	panel globals
	
		global	r4	0.2	
		global	r3	0.5	
		global	r2	1	
		global	r1	2	
		*ajg: i do this bc the #s in loop don't correspond to #s in panels

	///	set table globals
	
		global	post4	"\midrule"
		global	post3	"\midrule"
		global	post2	"\midrule"
		global	post1	""
		*ajg: i do this to avoid double line at the end of table
		
	foreach m in waz haz {
	
		eststo clear
		
		foreach v in `m'06 low_`m' very_low_`m' {
		
		///	all panels
	
			/// run analysis
			
				forval i=1/5 {
					qui eststo p1_`v'_w`i': areg wawc_`v'_c`i' treat ///
					bas_`m'06 $cov [aw=w] if bas_`m'06!=. & wawc_`m'06!=., ///
					absorb(stratum) cluster(digitcode) 
		
					qui su wawc_`v'_c`i' if treat==0 & bas_`m'06!=. & ///
					wawc_`m'06!=. 
					qui estadd scalar ymean=r(mean)
				}

			///	print table
			
				qui esttab p1_*_w5 using ${tables}estitt_`m'_win.tex, ///
				replace f refcat(treat "\emph{Winsorized at 0.1\%} (N=1538)", ///
				nolabel) nolines posthead(\midrule) booktabs label se ///
				drop(_cons $cov bas*) width(\hsize) b(%9.3f) nonotes ///
				star(* 0.1 ** 0.05 *** 0.01) stats(N ymean, fmt(%9.0f ///
				%9.3f %9.3f) labels("N (children)" "Control mean")) ///
				mlabels("\shortstack{${title_`m'}}" ///
				"\shortstack{${title_`m'2}}" "\shortstack{${title_`m'3}}") ///
				postfoot(\midrule)
			
				foreach r in 4 3 2 1 {					
					qui esttab  p1_*_w`r' ///
					using ${tables}estitt_`m'_win.tex, append f ///
					refcat(treat "\emph{Winsorized at ${r`r'}\%} (N=1538)", ///
					nolabel) nolines booktabs label se drop(_cons $cov ///
					bas*) width(\hsize) b(%9.3f) nonotes star(* 0.1 ** ///
					0.05 *** 0.01) stats(N ymean, fmt(%9.0f %9.3f %9.3f) ///
					labels("N (children)" "Control mean")) mlabels(none) ///
					nonumbers postfoot(${post`r'}) 	
				}
		}	
	}
					
///	table c1: raw proportion-correct scores at baseline
	
	/// load data 2
	
		use ${data}blanaldata2, clear
	
	///	print table
	
		balancetable treat per_math1 per_lang1 per_exec1 if in_bl==1 ///
		using ${tables}bas_propcorr.tex, replace booktabs leftctitle(" ") ///
		ctit("Control" "Treatment" "Difference" "N") observationscolumn ///
		noobs varlabels sd(par([ ])) fe(stratum) vce(cluster digitcode)
		
	
}
// ****************** GENERATE GRAPHS ****************** //

if $gengraphs==1 {

///	figure 1: benefit-cost sensitivity analysis

	set matsize 10000
	set seed 12345
	
	*****Set fixed parameter values
	set obs 1000000
	gen lfp=0.52
	gen wage=268
	gen te=0.11
	gen kids=14
	gen cohort_factor=1.33
	gen cost=74478
	gen te_haz=0
	local preferred=12.9
	gen g=0.14
	
	****Set random parameter distributions
	local days_lb=200
	local days_ub=250
	local growth_lb=0.03
	local growth_ub=0.07
	local r_lb=0.015
	local r_ub=0.045
	local mu_r=log(0.03)
	local sigma_r=0.694
	local p_lb=0.07
	local p_ub=0.19
	local p_haz_lb=0.20-(2*0.06)
	local p_haz_ub=0.20+(2*0.06)
	
	
	******Draw parameter values
	*Comment out third and fourth lines for uniform; comment out first line for trunc normal
	foreach t in r days growth p p_haz {
		gen `t' = ``t'_lb'+(``t'_ub'-``t'_lb')*uniform()
		local sigma=(``t'_ub'-``t'_lb')/4
		*gen `t'=``t'_lb'+((``t'_ub'-``t'_lb')/2)+`sigma'*invnormal(uniform())
		*drop if (`t'<``t'_lb') | (`t'>``t'_ub')
	}
	*Uncomment this for sensitivity analysis on discount rate
	*replace r=exp(`mu_r'+`sigma_r'*invnormal(uniform()))

	
	*****Compute PDV
	gen annual_earnings=lfp*wage*days
	foreach a of numlist 22/65 {
		gen pdv_`a'=(annual_earnings*((1+growth)^(`a'+4-5)))/((1+r)^(`a'-5))
	}
	egen pdv=rsum(pdv_*)
	gen benefit_kid=pdv*(p*(te - g*te_haz) + p_haz*te_haz)
	gen benefit=benefit_kid*kids*cohort_factor
	gen ratio=benefit/cost
	
	****Graph results
	sum ratio, detail
	local min=r(p5)
	local max=r(p95)
	keep if ratio<=90
	
		graph twoway kdensity ratio, ///
			 graphregion(color(white))  ///
			 ylabel(,nogrid) xline(`preferred', lpattern(dash) lcolor(black)) ///
			 xline(`min', lpattern(dash) lcolor(gray)) xline(`max', lpattern(dash) lcolor(gray)) ///
			 xlabel(0(10)90) ///
			xtitle("Benefit/cost ratio") ytitle("Density")
	
///	figure a3: children's age distribution by round of data collection

	/// load data 1
	
		use ${data}irtssanaldata, clear
	
	///	print graph
	
		tw (kdensity bas_child_age if in_bl==1 & (in_elawc==1|in_elhh==1) ///
		& std_math1!=. & bas_child_age!=., lcolor(blue) lpattern(solid)) ///
		(kdensity r5_child_age if in_bl==1 & in_elawc==1 & std_math2!=., ///
		lcolor(red) lpattern(longdash)) (kdensity hh_child_age [aw=w] if ///
		in_bl==1 & in_elhh==1 & std_math3!=. & (hh_sample==1| ///
		hh_sample==2), lcolor(green) lpattern(shortdash)), ///
		ytitle(Density) xtitle("Child age") legend(order(1 "Baseline" 2 ///
		"Center endline" 3 "Household endline") rows(1)) ///
		graphregion(fcolor(white) lcolor(white))

		gr export ${graphs}child_age.png, replace

///	figure a4: kling & liebman bounds

	forval x = 0(.01)0.3 {

	/// load data 1
	
		use ${data}irtssanaldata, clear
		
	///	#1: start w/kids from table 2
	
		/// sample indicator for awc
	
			g klsamp2=(fol_elawc!=. & treat!=. & female!=. & ///
			bas_child_age!=. & bas_waz06!=. & bas_haz06!=. & std_comb1!=.)
	
		///	sample indicator for hh
	
			g klsamp3=(fol_elhh!=. & treat!=. & female!=. & ///
			bas_child_age!=. & bas_waz06!=. & bas_haz06!=. & std_comb1!=.) & ///
			(hh_sample==1|hh_sample==2)

	///	#2-3: impute values for kids w/missing endline
	
		forval r=2/3 {
			foreach s in comb {				
				qui su stdf_`s'`r' if klsamp`r'==1 & treat==0
			
				qui replace stdf_`s'`r'=r(mean)+`x'*r(sd) if ///
				klsamp`r'==1 & treat==0 & stdf_`s'`r'==.
					
				qui su stdf_`s'`r' if klsamp`r'==1 & treat==1
				
				qui replace stdf_`s'`r'=r(mean)-`x'*r(sd) if ///
				klsamp`r'==1 & treat==1 & stdf_`s'`r'==.
			}
		}
		
	///	rerun analysis
	
		foreach s in comb {
			qui areg stdf_`s'2 treat std_`s'1 awc_std_`s'1 $cov ///
			if klsamp2==1, absorb(stratum) cluster(digitcode)
			
			*di _b[treat]
			di _se[treat]
		}	
	}
	
	///	plot graph
	
		tw (line comb_b x, lcolor(gray) yline(0, lcolor(black)) ///
		legend(order(1 "Composite score"))) (line comb_low x, ///
		lcolor(gray%50) lpattern(shortdash)) (line comb_upp x, ///
		lcolor(gray%50) lpattern(shortdash) graphregion(fcolor(white) ///
		lcolor(white)) xtit("Imputed values") ytit("Impact estimate"))
		
		graph export ${graphs}klbounds.png, replace
		
/// figure a5-a7: quantile effects

// *********** GRAPH QUANTILE EFFECTS - EDUCATION (BASELINE) ************** //

	
/// set title globals

	global	tit2	"Center assessments"
	global	tit3	"Household assessments"
	
///	set sampling restriction globals
	
	global 	sam2a	"in_elawc==1"
	global 	sam3a	"in_elhh==1"
 
	global 	sam2b	"in_elawc==1 & in_elhh==1"
	global 	sam3b	"in_elawc==1 & in_elhh==1"

///	set extension globals

	global	ext2	"awc"
	global	ext3	"hh"

	global 	exta	"full"
	global 	extb	"comm"
	
///	for each sampling restriction...	

	foreach s in a b {
		
///	foreach round...
	
	foreach r in 2 3 {
	
	/// load data 1
	
		use ${data}irtssanaldata, clear
		
		///	prepare dataset
	
			keep if in_bl==1 & ${sam`r'`s'}
		
		///	keep key vars
			
			keep digitcode childid stratum district treat std_comb*
		
		///	re-standardize scores w.r.t. sampling restriction
		
			su std_comb`r' if treat==0
			replace std_comb`r'=(std_comb`r'-r(mean))/r(sd)
	
		///	save as temp file
				
			tempfile temp_irt_data
			save `temp_irt_data'
		
		///	store iterations
			
			local it = 1000
			
		///	set seed
			
			set seed 12345
	
		///	bootstrap cis
			
			///	set range for plot
			
				qui egen kernel_range = fill(.01(.01)1)
		
			/// restrict kernel range to numbers below 1
			
				qui replace kernel_range = . if kernel_range>1
				
			///	convert kernel range to matrix if non-missing	
				
				mkmat kernel_range if kernel_range!=.
			
			/// input matrix for t-c diff
			
				matrix diff = kernel_range

			/// input matrix for matching criterion
				
				matrix x = kernel_range

			///	for each iteration...
			
				forvalues i = 1(1)`it'{
				
				/// ...call on data
					
					use `temp_irt_data', clear
	   
				///	...sample w/replacement
				
					bsample, strata(treat) cluster(digitcode)

				///	...rank scores by group
				
					bysort treat: egen rank=rank(std_comb1), unique
					
				///	...identify highest rank by group
				
					bysort treat: egen max_rank = max(rank)
					
				///	...identify percentile rank
				
					bysort treat: gen percentile = rank/max_rank

				///	...drop unecessary vars
				
					drop max_rank rank
		
				///	...set range for plot
				
					qui egen kernel_range = fill(.01(.01)1)
		
				/// ...restrict kernel range to numbers below 1
				
					qui replace kernel_range = . if kernel_range>1

				///	...regress endline score on percentile rankings
				
					lpoly std_comb`r' percentile if treat==0, ///
					gen(xcon_b2 y2_con_b) at (kernel_range) nograph
					lpoly std_comb`r' percentile if treat==1, ///
					gen(xtre_b2 y2_tre_b) at (kernel_range) nograph

				///	...convert to matrix if non-missing	
					
					mkmat y2_tre_b if y2_tre_b != . 
					mkmat y2_con_b if y2_con_b != . 
					matrix diff = diff, y2_tre_b - y2_con_b
				}
				
			///	transpose matrix
			
				matrix diff = diff'
				
			///	store columns from matrix as new variables
			
				svmat diff
			
			///	keep only variables from matrix
			
				keep diff*
				
			///	create matrix to store cis
			
				matrix conf_int = J(100, 2, 100)
				
			///	drop if original?
				*ajg: this does not seem to do anything
				
				qui drop if _n == 1
				
			///	sort each percentile and store 25th and 975th place in matrix
			
				/// for each iteration...
			
					forvalues i = 1(1)100{
					
					///	...sort differences
					
						 sort diff`i'
						 
					///	...store 25th and 975th place in matrix
					
						matrix conf_int[`i', 1] = diff`i'[25]
						matrix conf_int[`i', 2] = diff`i'[975]
					}
					
		///	graph data
		
			///	reload temp dataset
				
				use `temp_irt_data', clear
		
				///	rank WAZ scores by group
				
					bysort treat: egen rank=rank(std_comb1), unique
					
				///	identify highest rank by group
				
					bysort treat: egen max_rank = max(rank)
					
				///	identify percentile rank
				
					bysort treat: gen percentile = rank/max_rank

				/// drop unecessary vars
				
					drop max_rank rank
		
				///	set range for plot
				
					qui egen kernel_range = fill(.01(.01)1)
		
				/// restrict kernel range to numbers below 1
				
					qui replace kernel_range = . if kernel_range>1

				///	regress endline score on percentile rankings
				
					lpoly std_comb`r' percentile if treat==0, ///
					gen(xcon_b2 y2_con_b) at (kernel_range) nograph
					lpoly std_comb`r' percentile if treat==1, ///
					gen(xtre_b2 y2_tre_b) at (kernel_range) nograph
				
				///	g diff in coefficients
				
					gen diff = y2_tre_b - y2_con_b
		
				///	store columns from matrix as new variables
				
					svmat conf_int
	
				///	plot graph
				
					gr tw (line y2_con_b xcon_b2, lcolor(blue) ///
					lpattern("--.....") legend(lab(1 "Control"))) ////
					(line y2_tre_b xtre_b2, lcolor(red) ///
					lpattern(longdash) legend(lab(2 "Treatment"))) ///
					(line diff xcon_b2, lcolor(black) lpattern(solid) ///
					legend(lab(3 "Difference"))) (line conf_int1 xcon_b2, ///
					lcolor(black) lpattern(shortdash) ///
					legend(lab(4 "95% confidence band"))) ///
					(line conf_int2 xcon_b2, lcolor(black) ///
					lpattern(shortdash) legend(lab(5 ///
					"95% confidence band"))),yline(0, lcolor(gs10)) ///
					xtitle("Percentile of baseline composite score") ///
					ytitle("Endline composite score") ///
					legend(order(1 2 3 4)) title("${tit`r'}") ///
					graphregion(fcolor(white) lcolor(white)) ylab(-1(1)1)
					
					gr export ${graphs}qach_bas_${ext`r'}_${ext`s'}.png, ///
					replace	
	}
	}
	
// *********** GRAPH QUANTILE EFFECTS - EDUCATION (ENDLINE) ************** //


/// set title globals

	global	tit2	"Center assessments"
	global	tit3	"Household assessments"
	
///	set sampling restriction globals
	
	global 	sam2a	"in_elawc==1"
	global 	sam3a	"in_elhh==1"
 
	global 	sam2b	"in_elawc==1 & in_elhh==1"
	global 	sam3b	"in_elawc==1 & in_elhh==1"

///	set extension globals

	global	ext2	"awc"
	global	ext3	"hh"

	global 	exta	"full"
	global 	extb	"comm"
	
///	for each sampling restriction...	

	foreach s in a b {
		
///	foreach round...
	
	foreach r in 2 3 {
	
	/// load data 1
	
		use ${data}irtssanaldata, clear
		
		///	prepare dataset
	
			keep if in_bl==1 & ${sam`r'`s'}
			
		///	keep key vars
			
			keep digitcode childid stratum district treat std_comb`r'
		
		///	re-standardize scores w.r.t. sampling restriction
		
			su std_comb`r' if treat==0
			replace std_comb`r'=(std_comb`r'-r(mean))/r(sd)
	
		///	save as temp file
				
			tempfile temp_irt_data
			save `temp_irt_data'
		
		///	store iterations
			
			local it = 1000
			
		///	set seed
			
			set seed 12345
	
		///	bootstrap cis
			
			///	set range for plot
			
				qui egen kernel_range = fill(.01(.01)1)
		
			/// restrict kernel range to numbers below 1
			
				qui replace kernel_range = . if kernel_range>1
				
			///	convert kernel range to matrix if non-missing	
				
				mkmat kernel_range if kernel_range!=.
			
			/// input matrix for t-c diff
			
				matrix diff = kernel_range

			/// input matrix for matching criterion
				
				matrix x = kernel_range

			///	for each iteration...
			
				forvalues i = 1(1)`it'{
				
				/// ...call on data
					
					use `temp_irt_data', clear
	   
				///	...sample w/replacement
				
					bsample, strata(treat) cluster(digitcode)

				///	...rank scores by group
				
					bysort treat: egen rank=rank(std_comb`r'), unique
					
				///	...identify highest rank by group
				
					bysort treat: egen max_rank = max(rank)
					
				///	...identify percentile rank
				
					bysort treat: gen percentile = rank/max_rank

				///	...drop unecessary vars
				
					drop max_rank rank
		
				///	...set range for plot
				
					qui egen kernel_range = fill(.01(.01)1)
		
				/// ...restrict kernel range to numbers below 1
				
					qui replace kernel_range = . if kernel_range>1

				///	...regress endline score on percentile rankings
				
					lpoly std_comb`r' percentile if treat==0, ///
					gen(xcon_b2 y2_con_b) at (kernel_range) nograph
					lpoly std_comb`r' percentile if treat==1, ///
					gen(xtre_b2 y2_tre_b) at (kernel_range) nograph

				///	...convert to matrix if non-missing	
					
					mkmat y2_tre_b if y2_tre_b != . 
					mkmat y2_con_b if y2_con_b != . 
					matrix diff = diff, y2_tre_b - y2_con_b
				}
				
			///	transpose matrix
			
				matrix diff = diff'
				
			///	store columns from matrix as new variables
			
				svmat diff
			
			///	keep only variables from matrix
			
				keep diff*
				
			///	create matrix to store cis
			
				matrix conf_int = J(100, 2, 100)
				
			///	drop if original?
				*ajg: this does not seem to do anything
				
				qui drop if _n == 1
				
			///	sort each percentile and store 25th and 975th place in matrix
			
				/// for each iteration...
			
					forvalues i = 1(1)100{
					
					///	...sort differences
					
						 sort diff`i'
						 
					///	...store 25th and 975th place in matrix
					
						matrix conf_int[`i', 1] = diff`i'[25]
						matrix conf_int[`i', 2] = diff`i'[975]
					}
					
		///	graph data
		
			///	reload temp dataset
				
				use `temp_irt_data', clear
		
				///	rank WAZ scores by group
				
					bysort treat: egen rank=rank(std_comb`r'), unique
					
				///	identify highest rank by group
				
					bysort treat: egen max_rank = max(rank)
					
				///	identify percentile rank
				
					bysort treat: gen percentile = rank/max_rank

				/// drop unecessary vars
				
					drop max_rank rank
		
				///	set range for plot
				
					qui egen kernel_range = fill(.01(.01)1)
		
				/// restrict kernel range to numbers below 1
				
					qui replace kernel_range = . if kernel_range>1

				///	regress endline score on percentile rankings
				
					lpoly std_comb`r' percentile if treat==0, ///
					gen(xcon_b2 y2_con_b) at (kernel_range) nograph
					lpoly std_comb`r' percentile if treat==1, ///
					gen(xtre_b2 y2_tre_b) at (kernel_range) nograph
				
				///	g diff in coefficients
				
					gen diff = y2_tre_b - y2_con_b
		
				///	store columns from matrix as new variables
				
					svmat conf_int
	
				///	plot graph
				
					gr tw (line y2_con_b xcon_b2, lcolor(blue) ///
					lpattern("--.....") legend(lab(1 "Control"))) ////
					(line y2_tre_b xtre_b2, lcolor(red) ///
					lpattern(longdash) legend(lab(2 "Treatment"))) ///
					(line diff xcon_b2, lcolor(black) lpattern(solid) ///
					legend(lab(3 "Difference"))) (line conf_int1 xcon_b2, ///
					lcolor(black) lpattern(shortdash) ///
					legend(lab(4 "95% confidence band"))) ///
					(line conf_int2 xcon_b2, lcolor(black) ///
					lpattern(shortdash) legend(lab(5 ///
					"95% confidence band"))),yline(0, lcolor(gs10)) ///
					xtitle("Percentile of endline composite score") ///
					ytitle("Endline composite score") ///
					legend(order(1 2 3 4)) title("${tit`r'}") ///
					graphregion(fcolor(white) lcolor(white)) 
					
					gr export ${graphs}qach_${ext`r'}_${ext`s'}.png, ///
					replace	
	}
	}


// *********** GRAPH QUANTILE EFFECTS - NUTRITION (BASELINE) ************** //


/// set title globals

	global	tit_awc		"Center measurements"
	global	tit_hh		"Household measurements"
	
	global	tit_waz		"WAZ"
	global	tit_haz		"HAZ"
	
///	set sampling restriction globals
	
	global 	sam_awc_a	"in_elawc==1"
	global 	sam_hh_a	"in_elhh==1"

	global 	sam_awc_b	"in_elawc==1 & in_elhh==1"
	global 	sam_hh_b	"in_elawc==1 & in_elhh==1"

///	set extension globals

	global 	exta	"full"
	global 	extb	"comm"
	
///	foreach variable...

	foreach v in waz haz {
	
///	for each sampling restriction...	

	foreach s in a b {
	
///	foreach round...
	
	foreach r in awc hh {
			
	/// load dataset	
		
		use ${data}irtssanaldata, clear
		
		///	keep corresponding children
			
			keep if bas_`v'06!=. & ${sam_`r'_`s'}
				
		///	keep key vars
			
			keep digitcode childid stratum district treat bas_`v'06 `r'_`v'06
		
		///	save as temp file
				
			tempfile temp_irt_data
			save `temp_irt_data'
		
		///	store iterations
			
			local it = 1000
			
		///	set seed
			
			set seed 12345
		
		///	bootstrap cis
			
			///	set range for plot
			
				qui egen kernel_range = fill(.01(.01)1)
		
			/// restrict kernel range to numbers below 1
			
				qui replace kernel_range = . if kernel_range>1
				
			///	convert kernel range to matrix if non-missing	
				
				mkmat kernel_range if kernel_range !=.
			
			/// input matrix for t-c diff
			
				matrix diff = kernel_range

			/// input matrix for matching criterion
				
				matrix x = kernel_range

			///	for each iteration...
			
				forvalues i = 1(1)`it'{
				
				/// ...call on data
					
					use `temp_irt_data', clear
	   
				///	...sample w/replacement
				
					bsample, strata(treat) cluster(digitcode)

				///	...rank scores by group
				
					bysort treat: egen rank=rank(bas_`v'06), unique
					
				///	...identify highest rank by group
				
					bysort treat: egen max_rank = max(rank)
					
				///	...identify percentile rank
				
					bysort treat: gen percentile = rank/max_rank

				///	...drop unecessary vars
				
					drop max_rank rank
		
				///	...set range for plot
				
					qui egen kernel_range = fill(.01(.01)1)
		
				/// ...restrict kernel range to numbers below 1
				
					qui replace kernel_range = . if kernel_range>1

				///	...regress endline score on percentile rankings
				
					lpoly `r'_`v'06 percentile if treat==0, ///
					gen(xcon_b2 y2_con_b) at (kernel_range) nograph
					lpoly `r'_`v'06 percentile if treat==1, ///
					gen(xtre_b2 y2_tre_b) at (kernel_range) nograph

				///	...convert to matrix if non-missing	
					
					mkmat y2_tre_b if y2_tre_b != . 
					mkmat y2_con_b if y2_con_b != . 
					matrix diff = diff, y2_tre_b - y2_con_b
				}
				
			///	transpose matrix
			
				matrix diff = diff'
				
			///	store columns from matrix as new variables
			
				svmat diff
			
			///	keep only variables from matrix
			
				keep diff*
				
			///	create matrix to store cis
			
				matrix conf_int = J(100, 2, 100)
				
			///	drop if original?
				*ajg: this does not seem to do anything
				
				qui drop if _n == 1
				
			///	sort each percentile and store 25th and 975th place in matrix
			
				/// for each iteration...
			
					forvalues i = 1(1)100{
					
					///	...sort differences
					
						 sort diff`i'
						 
					///	...store 25th and 975th place in matrix
					
						matrix conf_int[`i', 1] = diff`i'[25]
						matrix conf_int[`i', 2] = diff`i'[975]
					}
					
		///	graph data
		
			///	reload temp dataset
				
				use `temp_irt_data', clear
		
				///	rank scores by group
				
					bysort treat: egen rank=rank(bas_`v'06), unique
					
				///	identify highest rank by group
				
					bysort treat: egen max_rank = max(rank)
					
				///	identify percentile rank
				
					bysort treat: gen percentile = rank/max_rank

				/// drop unecessary vars
				
					drop max_rank rank
		
				///	set range for plot
				
					qui egen kernel_range = fill(.01(.01)1)
		
				/// restrict kernel range to numbers below 1
				
					qui replace kernel_range = . if kernel_range>1

				///	regress endline score on percentile rankings
				
					lpoly `r'_`v'06 percentile if treat==0, ///
					gen(xcon_b2 y2_con_b) at (kernel_range) nograph
					lpoly `r'_`v'06 percentile if treat==1, ///
					gen(xtre_b2 y2_tre_b) at (kernel_range) nograph
				
				///	g diff in coefficients
				
					gen diff = y2_tre_b - y2_con_b
		
				///	store columns from matrix as new variables
				
					svmat conf_int
					
				///	plot graph
				
					gr tw (line y2_con_b xcon_b2, lcolor(blue) ///
					lpattern("--.....") legend(lab(1 "Control"))) ///
					(line y2_tre_b xtre_b2, lcolor(red) ///
					lpattern(longdash) legend(lab(2 "Treatment"))) ///
					(line diff xcon_b2, lcolor(black) lpattern(solid) ///
					legend(lab(3 "Difference"))) ///
					(line conf_int1 xcon_b2, lcolor(black) ///
					lpattern(shortdash) legend(lab(4 ///
					"95% confidence band"))) ///
					(line conf_int2 xcon_b2, lcolor(black) ///
					lpattern(shortdash) legend(lab(5 ///
					"95% confidence band"))) ,yline(0, lcolor(gs10)) ///
					xtitle("Percentile of endline ${tit_`v'} score") ///
					ytitle("Endline ${tit_`v'} score") ///
					legend(order(1 2 3 4)) ///
					title("${tit_`r'}") graphregion(fcolor(white) ///
					lcolor(white)) ylab(-4(1)1)
					
					gr export ${graphs}q`v'_bas_`r'_${ext`s'}.png, ///
					replace	
	}
	}
	}



// *********** GRAPH QUANTILE EFFECTS - NUTRITION (ENDLINE) ************** //


/// set title globals

	global	tit_awc		"Center measurements"
	global	tit_hh		"Household measurements"
	
	global	tit_waz		"WAZ"
	global	tit_haz		"HAZ"
	
///	set sampling restriction globals
	
	global 	sam_awc_a	"in_elawc==1"
	global 	sam_hh_a	"in_elhh==1"

	global 	sam_awc_b	"in_elawc==1 & in_elhh==1"
	global 	sam_hh_b	"in_elawc==1 & in_elhh==1"

///	set extension globals

	global 	exta	"full"
	global 	extb	"comm"
	
///	foreach variable...

	foreach v in waz haz {
	
///	for each sampling restriction...	

	foreach s in a b {
	
///	foreach round...
	
	foreach r in awc hh {
			
	/// load data 1
	
		use ${data}irtssanaldata, clear
		
		///	keep corresponding children
			
			keep if bas_`v'06!=. & ${sam_`r'_`s'}
				
		///	keep key vars
			
			keep digitcode childid stratum district treat bas_`v'06 `r'_`v'06
		
		///	save as temp file
				
			tempfile temp_irt_data
			save `temp_irt_data'
		
		///	store iterations
			
			local it = 1000
			
		///	set seed
			
			set seed 12345
		
		///	bootstrap cis
			
			///	set range for plot
			
				qui egen kernel_range = fill(.01(.01)1)
		
			/// restrict kernel range to numbers below 1
			
				qui replace kernel_range = . if kernel_range>1
				
			///	convert kernel range to matrix if non-missing	
				
				mkmat kernel_range if kernel_range !=.
			
			/// input matrix for t-c diff
			
				matrix diff = kernel_range

			/// input matrix for matching criterion
				
				matrix x = kernel_range

			///	for each iteration...
			
				forvalues i = 1(1)`it'{
				
				/// ...call on data
					
					use `temp_irt_data', clear
	   
				///	...sample w/replacement
				
					bsample, strata(treat) cluster(digitcode)

				///	...rank scores by group
				
					bysort treat: egen rank=rank(bas_`v'06), unique
					
				///	...identify highest rank by group
				
					bysort treat: egen max_rank = max(rank)
					
				///	...identify percentile rank
				
					bysort treat: gen percentile = rank/max_rank

				///	...drop unecessary vars
				
					drop max_rank rank
		
				///	...set range for plot
				
					qui egen kernel_range = fill(.01(.01)1)
		
				/// ...restrict kernel range to numbers below 1
				
					qui replace kernel_range = . if kernel_range>1

				///	...regress endline score on percentile rankings
				
					lpoly `r'_`v'06 percentile if treat==0, ///
					gen(xcon_b2 y2_con_b) at (kernel_range) nograph
					lpoly `r'_`v'06 percentile if treat==1, ///
					gen(xtre_b2 y2_tre_b) at (kernel_range) nograph

				///	...convert to matrix if non-missing	
					
					mkmat y2_tre_b if y2_tre_b != . 
					mkmat y2_con_b if y2_con_b != . 
					matrix diff = diff, y2_tre_b - y2_con_b
				}
				
			///	transpose matrix
			
				matrix diff = diff'
				
			///	store columns from matrix as new variables
			
				svmat diff
			
			///	keep only variables from matrix
			
				keep diff*
				
			///	create matrix to store cis
			
				matrix conf_int = J(100, 2, 100)
				
			///	drop if original?
				*ajg: this does not seem to do anything
				
				qui drop if _n == 1
				
			///	sort each percentile and store 25th and 975th place in matrix
			
				/// for each iteration...
			
					forvalues i = 1(1)100{
					
					///	...sort differences
					
						 sort diff`i'
						 
					///	...store 25th and 975th place in matrix
					
						matrix conf_int[`i', 1] = diff`i'[25]
						matrix conf_int[`i', 2] = diff`i'[975]
					}
					
		///	graph data
		
			///	reload temp dataset
				
				use `temp_irt_data', clear
		
				///	rank scores by group
				
					bysort treat: egen rank=rank(bas_`v'06), unique
					
				///	identify highest rank by group
				
					bysort treat: egen max_rank = max(rank)
					
				///	identify percentile rank
				
					bysort treat: gen percentile = rank/max_rank

				/// drop unecessary vars
				
					drop max_rank rank
		
				///	set range for plot
				
					qui egen kernel_range = fill(.01(.01)1)
		
				/// restrict kernel range to numbers below 1
				
					qui replace kernel_range = . if kernel_range>1

				///	regress endline score on percentile rankings
				
					lpoly `r'_`v'06 percentile if treat==0, ///
					gen(xcon_b2 y2_con_b) at (kernel_range) nograph
					lpoly `r'_`v'06 percentile if treat==1, ///
					gen(xtre_b2 y2_tre_b) at (kernel_range) nograph
				
				///	g diff in coefficients
				
					gen diff = y2_tre_b - y2_con_b
		
				///	store columns from matrix as new variables
				
					svmat conf_int
					
				///	plot graph
				
					gr tw (line y2_con_b xcon_b2, lcolor(blue) ///
					lpattern("--.....") legend(lab(1 "Control"))) ///
					(line y2_tre_b xtre_b2, lcolor(red) ///
					lpattern(longdash) legend(lab(2 "Treatment"))) ///
					(line diff xcon_b2, lcolor(black) lpattern(solid) ///
					legend(lab(3 "Difference"))) ///
					(line conf_int1 xcon_b2, lcolor(black) ///
					lpattern(shortdash) legend(lab(4 ///
					"95% confidence band"))) ///
					(line conf_int2 xcon_b2, lcolor(black) ///
					lpattern(shortdash) legend(lab(5 ///
					"95% confidence band"))) ,yline(0, lcolor(gs10)) ///
					xtitle("Percentile of endline ${tit_`v'} score") ///
					ytitle("Endline ${tit_`v'} score") ///
					legend(order(1 2 3 4)) ///
					title("${tit_`r'}") graphregion(fcolor(white) ///
					lcolor(white)) ylab(-4(1)1)
					
					gr export ${graphs}q`v'_bas_`r'_${ext`s'}.png, ///
					replace	
	}
	}
	}
		
///	figure a8: benefit-cost sensitivity analysis (log-normal discount rate)

	set matsize 10000
	set seed 12345
	
	*****Set fixed parameter values
	set obs 1000000
	gen lfp=0.52
	gen wage=268
	gen te=0.11
	gen kids=14
	gen cohort_factor=1.33
	gen cost=74478
	gen te_haz=0
	local preferred=12.9
	gen g=0.14

	
	****Set random parameter distributions
	local days_lb=200
	local days_ub=250
	local growth_lb=0.03
	local growth_ub=0.07
	local r_lb=0.015
	local r_ub=0.045
	local mu_r=log(0.03)
	local sigma_r=0.694
	local p_lb=0.07
	local p_ub=0.19
	local p_haz_lb=0.20-(2*0.06)
	local p_haz_ub=0.20+(2*0.06)
	
	
	/*Pick model 
		*truncnorm_params=0 and lognorm_r=0 for fig 1a
		*truncnorm_params=1 and lognorm_r=0 for fig 1b
		*truncnorm_params=0 and lognorm_r=1 for fig a8a
		*truncnorm_params=1 and lognorm_r=1 for figa8b
	*/
	local truncnorm_params=0
	local lognorm_r=0
	
	
	******Draw parameter values
	*Comment out third and fourth lines for uniform; comment out first line for trunc normal
	foreach t in r days growth p p_haz {
		if `truncnorm_params'==0 {
			gen `t' = ``t'_lb'+(``t'_ub'-``t'_lb')*uniform()
		}
		if `truncnorm_params'==1 {
			local sigma=(``t'_ub'-``t'_lb')/4
			gen `t'=``t'_lb'+((``t'_ub'-``t'_lb')/2)+`sigma'*invnormal(uniform())
			drop if (`t'<``t'_lb') | (`t'>``t'_ub')
		}
	}
	if `lognorm_r'==1 {
		replace r=exp(`mu_r'+`sigma_r'*invnormal(uniform()))
	}
	
	*****Compute PDV
	gen annual_earnings=lfp*wage*days
	foreach a of numlist 22/65 {
		gen pdv_`a'=(annual_earnings*((1+growth)^(`a'+4-5)))/((1+r)^(`a'-5))
	}
	egen pdv=rsum(pdv_*)
	gen benefit_kid=pdv*(p*(te - g*te_haz) + p_haz*te_haz)
	gen benefit=benefit_kid*kids*cohort_factor
	gen ratio=benefit/cost
	
	****Graph results
	sum ratio, detail
	local min=r(p5)
	local max=r(p95)
	
	*For graphical presentation, eliminate ratios > 90 for sensitivity analysis
	if `lognorm_r'==1 {
		keep if ratio<=90
	}
	
		graph twoway kdensity ratio, ///
			 graphregion(color(white))  ///
			 ylabel(,nogrid) xline(`preferred', lpattern(dash) lcolor(black)) ///
			 xline(`min', lpattern(dash) lcolor(gray)) xline(`max', lpattern(dash) lcolor(gray)) ///
			 xlabel(0(10)90) ///
			xtitle("Benefit/cost ratio") ytitle("Density")

///	figure c1: raw proportion-correct scores in assessments by round

	/// set title globals
	
		global tit_math	"Math"
		global tit_lang	"Language"
		global tit_exec	"Executive function"
		global tit_comb	"Composite"

	/// load data 1
	
		use ${data}irtssanaldata, clear
		
	/// print graph
	
		foreach s in $subjects {
			tw (hist per_`s'1 if stdf_`s'2!=. & std_`s'1!=. & ///
			aww_pass10!=. & m_pass10!=. & aww_exper!=. & m_exper!=., ///
			percent width(.04) color(navy%50)) ///
			(hist per_`s'2 if stdf_`s'2!=. & std_`s'1!=. & ///
			aww_pass10!=. & m_pass10!=. & aww_exper!=. & m_exper!=., ///
			percent width(.04) color(green%50)) ///
			(hist per_`s'3 if stdf_`s'3!=. & std_`s'1!=. & ///
			aww_pass10!=. & m_pass10!=. & aww_exper!=. & m_exper!=., ///
			percent width(.04) color(maroon%50)), ylab(0(20)80) ///
			tit(${tit_`s'}) legend(order(1 "Baseline" 2 "Center endline" ///
			3 "Household endline") rows(1)) ///
			graphregion(fcolor(white) lcolor(white))
			gr save ${graphs}per_`s'.gph, replace
		}
		
		grc1leg ${graphs}per_math.gph ${graphs}per_lang.gph ///
		${graphs}per_exec.gph ${graphs}per_comb.gph, ///
		leg(${graphs}per_comb.gph) ///
		graphregion(fcolor(white) lcolor(white)) 
		gr export ${graphs}dists_per.png, replace
	
}
